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INTRODUCTION 

In problems of impact such as those experienced by a vehicle landing on 
some planet, entering into some ocean or subjected to some high rate distur- 
bances of the external pressure, the compressive stress wave emanating from 
the surface of contact could be quite severe. This may cause severe damage 
or failure to the vehicle or to the various components within the vehicle. 
Therefore, it becomes necessary to mitigate the loads felt by the payload 
so that the maximum allowable stress is not exceeded. The purpose of this 
effort is to study both analytically and experimentally the attenuation of 
the stress waves during passage through single and multilayer structures. 

One method of attenuating the magnitude of the incident compressive 
stress is to incorporate into the system a material or multilayered material 
which may be permitted to deform plastically, thereby mitigating the load. 

The problem has been treated thus far with reasonable success by empirical 
means, but the intelligent design of such a shock mitigating system requires 
the understanding and the application of the theory of elastic and plastic 
wave propagation as well as' extensive knowledge of the material properties 
of the mitigating device. With the support from the National Aeronautical 
and Space Administration to which this final report is submitted, several 
basic research studies, both analytically and experimentally, have been 
carried out with conclusive results. 

The investigation included studies on elastic and plastic stress wave 
propagation in the composites and those on shock mitigating material character- 
istics such as dynamic stress-strain relations and energy absorption properties. 
Studies were performed on the behavior of composites under dynamic loading. 

A summary of the studies is presented in the following: 
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I. Stress wave propagation in composites: 

A. Propagation of elastic waves in composites. 

B. Propagation of elastic waves through areas of discontinuity. 

C. Propagation of radial elastic plastic waves in a circular disk. 

D. Elastic and Plastic wave cancel lation. 

II. Shock mitigating material characteristics: 

A. Energy absorption characteristics of several shock mitigating materials. 

B. Dynamic stress-strain relations of several shock mitigating materials. 

III. Dynamic response of composites: 

A. The longitudinal vibration of axially symmetric composite bodies. 

B. Shock mitigation models. 
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I. STRESS WAVE PROPAGATION IN COMPOSITES 
A. Propagation of Elastic Waves in Composites 

The propagation of an initally sharp cylindrical pressure pulse through 
a linear elastic fiber reinforced composite medium is analysed, both experi- 
mentally and analytically. In the experiment, tests were performed on plates 
with single and multiple circular inclusions embedded in a matrix of lower 
characteristic impedance. Sharp compression pulses were generated at an 
edge of the plate by two methods: the detonation of a charge of lead azide 

and the impingement of a sharp knife-edged ram impacted at the opposite end 
by a short projectile. Strain gages were mounted on various positions of 
the plate to determine the attenuation of the transient stress in the fiber 
reinforced composite. The qualitative analytical treatment is based on the 
methods of propagating stress discontinuities. Computer programs were written 
to numerically determine the changes in the shape of the leading wave front 
and the stresses immediately behind it. Experimental results for the attenua- 
tion of stress wave on steel -a 1 uminum and steel -brass fiber-matrix composites 
compared very well with the computed analytical results when the applied 
pressure is generated by small explosive charges. The results did not compare 
well when the applied pressure is generated by projectile impact. 

Results of this portion of the investigation were presented and published 
at the AIAA/ASME/SAE 13th Structures, Structural Dynamics, and Materials 
Conference in San Antonio, Texas, April 1972. The paper is entitled, "Experi- 
mental and Theoretical Investigation of Stress Wave Attenuation in Fiber 
Reinforced Composites." A copy of the paper is included in the Appendixes. 

-B. -Propagation of Elastic Waves Through Areas of Discontinuity 

Elastic Stress Waves generated in structure consisting of cylindrical 
rods with grooves, holes, threads, and screws by the application of pressure 
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pulses of various durations are investigated both experimentally and analytically. 
In the experiment, the stress pulse was introduced by using an air gun to 
impact aluminum projectiles on various configured aluminum rods. A schematic 
diagram of the experimental set-up is shown in Fig. 1. The 3/4 inch aluminum 
rods were machined to the configurations shown in Fig. 2 and 3* Strain gages 
were mounted on the rod in diametrically opposed pairs at locations shown in 
Figs. 2 and 3. Each pair of gages was connected to a full Wheatstone bridge 
using two dummy gages so that bending effects would be cancelled. The effect 
of grooves, holes, threads and screws on elastic waves in cylindrical rods 
were obtained from the various parameters from the experimental tests. A 
typical result is given in Fig. 4 for tests on a cylindrical rod with a 2 inch 
deep hole at the opposite end of impact. The maximum stress measured from 
the strain gage is obtained for various hole diameter. This is repeated for 
plain hole, hole with thread and hole with thread and screw. Another result 
is given in Fig. 5 for tests on a cylindrical rod with a 1 inch deep hole at 
the opposite end of impact. This and other results indicated that the threads 
and grooves had very little effect on elastic wave propagation. The results 
are very realistic and very predictable. 

In the theoretical analysis, both a one-dimensional and a two-dimensional 
nimterical method for the solution of elastic waves in various configured rods 
are utilized. For the one dimensional analysis the method of characteristics 
is employed to analyze the propagation of stress waves in the rods. A 
numerical integration of the characteristic equation in each section, in 
conjunction with appropriate procedures to ensure that the conditions on the 
external boundaries as well 'as the continuity conditions at the intersections 
are satisfied, yields the stress and the particle velocity of any arbitrary 
location. For the two-dimensional analysis, the exact equations of motion 
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FIGURE I SCHEMATIC DIAGRAM OF INSTRUMENTAL SETTING 
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. FIGURE 4 ..CYLINDRICAL ROD WITH 2" DEEP HOLE 
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governing elastic, axisymmetric wave propagation in a cylindrical rod are 
approximated by a first-order finite-difference scheme. This difference 
scheme is based on a displacement rather than a velocity formulation, 
thereby making it unnecessary to explicitly introduce an artificial viscosity 
term into the finite-difference equations. Numerical results are obtained 
using the resulting difference equations. in conjunction with the boundary 
and initial conditions. The numerical results compared very well with the 
experimental data. From a design standpoint, it is very desirable to be 
able to predict the magnitude of the transmitted and reflected protions of 
the waves . 

Results of this protion of the investigation constituted the basic 
material for a Master's thesis and a paper submitted for presentation in 
the S.E.S.A. Winter Annual Meeting in Seattle, Washington, Oct. 1972. A 
copy of the thesis is included in the Appendixes. 

C. Propagation of Radial Elastic-Plastic V/aves in a Circular Disk. 

The prupose of this work was to investigate the problem of radial 
elastic-plastic stress wave propagation in a thin, finite, work hardening 
disk using the method of characteristics and the method of successive elastic 
approximations. In this method the characteristic equations, derived for 
the elastic case, are applied to the plastic regime by allowing the modulus 
of elasticity E and the wave speed C appearing in the equations to assume 
values corresponding to the tangent modulus of the plastic stress-strain curve. 
The value of Poisson 1 Ratio, is taken as 1/2 in the plastic range. 

Three characteristic equations were derived from the governing conserva- 
tion and elastic constituent equations as outlined in a paper by Chou and 
Koenig. They were used subsequently, along with the Prandt 1 -Reuss incremental 
plasticity theory and the method of successive elastic solutions, to calculate 



the problem variables at the nodes of the characteristic network which was 
generated in a step-by-step fashion as the solution progressed. 

A modest experimental program consisting of two experimental techniques 
was conducted to check the accuracy of the proposed analytical method. A 
circular disk was used with a hole at the center to accommodate a pressure 
pulse applied either hydraulically or explosively to the inner surface. 

The hydraulic method produced poor results and was di sea rded in favor of 
the explosive technique. 

The disk material was 1100-0 aluminum which exhibits a 0.2 percent 
offset yield strength of 5000 psi and a proportional limit of 2500 psi. 
Strain-rate effects were not included in the theory. 

Two 350 ohm Micro-Measurement strain gages were used in addition to a 
TEKTRONIX Model 551 dual beam osci 1 loscope to obtain radial and circum- 
frential strain time histories at a radial distance of 0-75 inches on a 
k inch 0-D. specimen. 

A Piezotronics Model 109A ballistics pressure transducer was used 
with one channel of a second Model 551 TEKTRONIX dual beam oscilloscope 
to record the applied pressure pulse. A computer program incorporating 
the analytical procedures described previously was used to check the strain 
histories measured on the specimens. In addition, the theoretical solutions 
for a ramp followed by a steady state pressure were run as approximate 
checks against the known, closed form, static elastic and elastic-plastic 
solutions for a thick-walled cylinder. 

Reasonable agreement between the theoretical and experimental results 
supports the hypothesis that this method can be used to solve problems of 
wave propagation in systems subjected to combined states of stress. 
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Results of this portion of the investigation constituted the basic 
material for a Ph.D. Thesis. A copy of the thesis is included in the 
Appendixes. 

D. Elastic and Plastic Wave Cancellation 

In order to understand this concept, it should be recalled that 
after a material has been loaded to a point where plastic deformation 
occurs., while the load is being released, the displacement will obey the 
Hooke's law instead of retracting its original path. Since stress waves 
will propagate at a speed directly proportional to the square root of the 
slope of the stress-strain curve, it may be seen that the transmission of 
stress waves during the unloading will take;palce much faster than the 

V . 

transmission of the loading stresses above the yield point. Therefore, this 
method of mitigating the applied impulse is to make use of the overtaking 
and cancellation of high-intensity plastic stress waves by the faster-traveling 
elastic unloading stress waves. 

A study was made to further develop the one-dimensional stress wave 
theory and confirm the method by which stress attenuation is attained. A 
one-dimensional elastic-plastic code was developed for the calculations of 
the stress deformation and particle velocity history. This code is applicable 
to problems of one-dimensional longitudinal wave propagation in rods and 
radial stress-wave propagation in circular disks. The method of characteristics 
was used and the characteristic equations were derived from the governing 
equilibrium and elastic constituent equations. These equations were used 
subsequently along with an incremental plasticity theory embodying the method 
of successive elastic approximations, to calculate the problem variables 
at the nodes of the characteristic network which was generated in a step-by- 
step fashion as the solution progresses. 
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A series of experimental tests were performed in order to verify 
the analytical results. One such experiment utilizes an air gun to 
accelerate a steel rod and impact it against an instrumented test rod protected 
by various mitigating materials or by various multilayered mitigating materials. 
A stress wave was developed in the test rod which was measured with strain 
gages. The reduced stress obtained was analyzed as a function of the applied 
force pulse, the material properties of the mitigators., the length of the 
mitigators, and the number and combination of the mitigators for optimum 
mitigation. The results, as shown in Figures 6 and 7 indicated that the 
polymer group, such as ethocel, teflon, lexan, and polyurethane, is very 
efficient for shock mitigation from both a weight and volume standpoint. 

Then to verify the concept of wave cancellation, the theoretical 
determination of the stress and the strain distributions in a test where 
an air gun was used to accelerate a lexan rod and impact it against a second 
stationary lexan rod placed in front of an elastic measuring bar was performed. 
The stress and strain distributions obtained in such a dynamic compression 

f 

test compared very satisfactorily with the stress and the strain distributions 
deduced from the theoretical analysis. Figures 8 through 10 give the comparison 
of the strain measurements from the strain gages to the theoretically calculated 
values for various lengths of mitigators with lexan projectiles of l/2, 1, 
and 2 inches in length, respectively. Figure 11 gives the comparison of 
the strain measurements from the strain gages to the theoretically calculated 
values with a steel projectile 1 inch in length. From these results, it 
could be recommended that for design purposes, a mitigator 1 ength-to-proj ect i 1 e 
length ratio of 0.5 seems to offer the optimum shock mitigation for the 


materials tested. 
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II. SHOCK MITIGATING MATERIAL CHARACTERISTICS 

A. Energy Absorption Characteristics of Several Shock Mitigating Materials 

The scope of this investigation is limited to the determination of the 
absorption characteristic with respect to the thickness of the shock mitigating 
materials under low impact velocity conditions. The materials tested included 
styrofoam, polyethylene, high and low density polyurethene and Balsa wood. 

Specimens 3 /4 inches in diameter were placed coaxially at the muzzle 
of the 3/4 inch gun with one end in contact with an instrumented rod, (Fig. 12). 

The impact was initiated by a low velocity bullet crushing against the specimen. 

The peak stress measured with strain gages in the instrumented rod was compared 
with that resulted in direct impact of the bullet to the instrumented rod. 

Similar experiments were also conducted by means of a 6 inch dia. drop weight 
tester, (Fig. 13). In this case, the specimens were crushed between a falling 
cylindrical weight and a stationary instrumented platform. 

Experimental results indicated that for styrofoam, ployethylene and Balsa 
wood, the shock absorption capability generally increases with thickness. On 
the other hand, polyurethene was found to be a very poor shock absorbing 
material. The high density polyurethene retained virtually the same thickness 
under impact. The low density polyurethene was compressed easily but bounced 
back to its original thickness while offering little shock mitigating capa- 
bility up to as thick as 4 inches. 

Shock mitigating characteristics for styrofoam, polyethene and Balsa 
wood are presented in Figures 14 through 17* 

B. Dynamic Stress-Strain Relations of Several Shock Mitigating Materials 

Lightweight foam materials are used abundantly in contemporary design 
for shock isolation. It is well known that the dynamic properties of many 
materials differ from the static properties. Tests were performed to obtain 
dynamic stress-strain data for a representative foam material -rigid polyurethane foam 




Figure 12 Air Gun Impact Test Apparatus. 
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Figure 13 6" Drop Weight Tester. 


19 



uares tu the Inch 










20 



i Thickness ' in I Inches' 










































21 























































23 


Three densities of the foam (nominal 3*5j 7-0 and 14.5 lb/ft^) were 

tested in compression and tension over a wide range of strain rates. The 

-2 4 

rates were varied between 10 and 10 in/in/sec for the compression tests 
"*2 

and between 10 and 10 in/in/sec for the tension tests. The foam was 
tested parallel to the direction of rise and perpendicular to the direction 
of rise, to check for anisotropy. 

Stress-strain relations .were thus obtained as functions of strain rate 
for each of the densities and rise directions. Yield strengths and elastic 
moduli were then determined as functions of strain rate from the compressive 
stress-strain data. No definite yield points were observed in tension. 

Therefore, ultimate strengths and elastic moduli were obtained as functions 

ZTN 

of strain rage from the tensile stress-strain data. 

W 

Two test apparatus were required to obtain data over the wide range of 

strain rates. A Plastechon Model 581 Universal testing machine was used to 

-2 

produce the rates of strain between 10 and 10 in/in/sec. A split Hopkinson 
pressure bar was used to produce the rates of strain between 10^ and 10^ in/in/sec 
in compression. 

One inch diameter cylindrical specimens having lengths of 0.5, 1 and 2 
inches and axes parallel and perpendicular to the foam rise direction were 
tested in compression using the Model 581 machine. The 0.4 inch diameter 
cylindrical specimens having lengths of 0.2, 0.4, and 0.7 inches were 
tested in compression on the split Hopkinson bar. Specimens having 0.375 
inches square cross-sections and 1.625 inch gauge lengths were tested in 
tension on the Model 581 machine. 

Conclusions are drawn that nominal 3-5> 7*0 and 14.5 lb/ft J rigid 
polyurethane foams are strain rate sensitive in compression. The yield 

-2 

strength increases consistently over the range of strain rates between 10 



and 10 in/in/sec. The ultimate tensile strength and tensile elastic modulus 
are also not rate dependent between 10 and 10 in/in/sec. Any anisotropic 
behavior is small and obscurred by the scatter in the data. An increase in 
foam density results in an increase in both the yield strength, ultimate 
strength and elastic modulus. 
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III. DYNAMIC RESPONSE OF COMPOSITES 

A. The Longitudinal Vibration of Axially Symmetric Composite Bodies. 

The natural frequencies and corresponding mode shapes for composite, 
axisymmetric bodies of varying area undergoing longitudinal vibrations are 
determined both experimentally and theoretically. Composite rods, cones and 
conical shells of Lucite, Nylon, Polycarbonate and Polyethylene are investi- 
gated. .Strain measurements are taken from strain gages mounted on the surface 
of the models which are excited at varying frequencies by a shaker table. 

Lumped parameter approximations of the continuous bodies are solved numerically 
on the digital computer for both fixed-free and free-free boundaries in order 
to determine the natural frequencies, mode shapes and relative strains. The 
method of character isitcs is utilized in a wave propagation approach to obtain 
the response of a point in the model to a random input. The random input 
selected is white noise. The natural frequencies are obtained by performing 
a Fast Fourier Transform Analysis on the response. 

Results of this portion of the investigation were published in the 
Shock and Vibration Bulletin, Bulletin hi, Part 5, pp. 227-23*+, January 1972. 

The paper is entitled "Longitudinal Vibration of Composite Bodies of Varying 
Area". A copy of the paper is included in the Appendixes. 

B. Shock Mitigation Model 

Many mater ia 1 s, e.g., rigid polyurethane foam, exhibit elastic, perfectly 
plastic stress-strain curves until compaction occurs at large strains. Stresses 
greater than the yield stress cannot be transmitted by such materials until 
compaction occurs. An investigation has been initiated using a lumped parameter 
spring mass analysis (for the constant area case) and a multiple spring mass 
analysis (for the changing area case), see Figure 18. 
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The mitigator is divided into a series of masses connected by springs 
and dashpots, figure 18. The force-deflection characterisitcs of the springs 
are related to the stress-strain relation of the mitigating material and the 
cross-sectional area of the corresponding mass.- A rigid, perfectly plastic 
or linearly work-hardening, locking stress-strain curve is assumed for the 
material, Figure 19* No deflection recovery is allowed upon unloading. The 
dashpot is considered a locking device and can be used to incorporate strain 
rate effects . 

The following experiments were performed to verify the theoretical analysis. 
The tests consist of impacting cylindrical and conical mitigators with cylindrical 
projectiles and measuring the transmitted force, Figure 20. The stress and 
strain distributions obtained in such dynamic compression tests will be compared 
with those calculated from the theoretical analysis. Two of the test results 
of the transmitted strain measurements are shown in Figures 21 and 22 for 
cylindrical and conical mitigators, respectively. It was observed that the 
transmitted strain did not remain at a constant level as expected for the 
constant area case. This drop in strain could be due to the strain-rate 
effect, since the rate of strain decreases as the mitigator is being crushed, 
thereby lowering the yield strength of the mitigator. 

The projectile energy was varied by varying both the projectile mass 
and velocity. The energy absorbtion capability of the foam was varied by 
varying the length of the mitigator. In the future, the energy absorption 
capability will be varied. by varying the foam density. 
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ABSTRACT 


Title of Thesis: THE PROPAGATION OF ELASTIC STRESS WAVES IN 

DISCONTINUOUS STRUCTURES 
Supachai Limpisvasti, Master of Science, 1972 
Thesis directed by: Drs. C. Y. Tsui and J. C. S. Yang 

The propagation of elastic stress waves of various pulse durations in 
cylindrical rods with grooves, holes, threaded holes and threaded holes 
with screws were investigated both experimentally and analytically. 

In the experimental tests, the input stress pulse was generated in 
the configured aluminum rod by the impact of a projectile issuing from 
an air gun. Dynamic strain gages were used as the sensing device. 

Analytical results were obtained by the method of characteristics 
for one dimensional wave propagation. 

• The effects of various sizes of grooves, holes, etc. on elastic 

• • •'.*•* j ... 

stress wave propagation were examined and the decay of stress in regions 
close to the section of discontinuity was also studied. 
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INTRODUCTION 


Recently, increasing interest has developed in the transient respors 
of structural elements involving the propagation of elastic stress waves. 
Among the different methods used in these investigations, a numerical 
method based on the integration of the characteristic equations along 
characteristic lines has proved to be simple and accurate. Chou with 
Mortimer (1) used this method to find the solution of the propagation 
of discontinuities in one dimensional plane elastic wave and then with 
Koeing (2) to find solutions of a step input stress for one dimensional 
cylindrical and spherical elastic waves. Also, Meng and McNiven ( 3 ) 
-analysed the transient excitation of an elastic rod by this method. 

The method of characteristics is not restricted to single layer, 
homogeneous, constant area, elastic structures. Tang and Hasset (U) 
investigated stress waves in multilayered axisymmetrical bodies of 
varying area. Mao and Rader (5) investigated the longitudinal stress 
propagation in non-uniform elastic and viscoelastic bars. 

The purpose of this thesis is to investigate both experimentally 
and analytically the propagation of elastic stress waves due to an impact 
loading in discontinuous structures such as rods vrith grooves, holes, 
threaded holes and threaded holes with screws. 

For the analytical approach, the method of characteristics for 

t 

one dimensional wave propagation is employed. 

The distinguishing feature of impulsive loading arises from the 
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fact that the action of a suddenly applied load is not simultaneously 
felt at all parts of the body. Parts of the body remote from the point 
of application of the load remain undisturbed until the stress wave 
in reaction to the load is propagated to those parts. The finite 
velocity of propagation of such waves depends upon' inertial effects 
and upon mechanical properties of the material comprising the body. 

The governing equation of the propagation of stress waves can be 
developed by means of either the Eulerian or the Lagrangian coordinate 
system. The Lagrangian coordinate system, however, exhibits one 
significant advantage over the Eulerian coordinate system in that the 
former provides a simpler expression for the law of conservation of mass. 
Therefore, in this thesis, the governing equation will be expressed in 
the Lagrangian coordinate system. 

In the experimental investigations, strain gages were used as sensing 
devices to measure the strain caused by the propagation of the elastic 
stress wave in the specimen. Strain gages have been used by many investi- 
gators in such studies, giving very reliable results especially in the 
elastic range. Bickle (6) published the results of the use of strain 
gages for the measurement of propagating strain waves and applied them 
to mary cases. Valather and Baker (7) used strain gages in investigating 
-wave propagation and obtained good results. The input loading pulse is 
obtained by impacting various aluminum projectiles against the specimen 
with the aid of an air gun. The results will be analyzed and compared 
with analytical solutions. 
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CHAPTER I 

THEORETICAL SOLUTION 

In the theoretical analysis of stress wave propagation in the elastic 
region, three methods of solution are generally used, namely, the 
Laplace, transformation method, the mode superposition method and the 
method of characteristics. In this thesis, the method of characteristics 
will be employed because the mathematics involved in the solution is 
simple and because the method can be applied very simply to any arbitrary 
input wave form. 

The problem is easily described in the Lagrangian coordinate 
•system with a single spatial variable, a, referred to the initial 
position of each layer of particles. In this system, the co-ordinate 
moves along with the particles which results in certain mathematical 
simplifications. In addition, since many strain measuring devices are 
attached to a set of particles and move along with these particles, the 
strain measurements are directly related to those predicted from this 



Inertial Reference Line total 


A* Method of Characteristics 


The particle velocity is defined as: 
u - f^x 


w\ 


where x is the Eulerian coordinate, t is time and a is the Lagrangian 
coordinate. Since the engineering strain is caused by the change in x 
due to a change in a, it can be defined as: 

(!?)t da - da 


da 


or 


(&) * = i + * w 

For an infinitesimal alice of an object, the mass, M, will be the 
same for the strained and the unstrained conditions: 


M - D da * 


A 


'©~o 




P-*\ _A A n' 
Uajt -S-fit 


or 


$ 


O^o 


/* 

- 1 +£ 


( 2 ) 


Also, the rate of mass, m, entering the wave front will be equal to 
that leaving the wave front: 

m = C/3 A - (C - du)(/OA + d(/oA) ) 

where C is the velocity of proapgation of the wave front with respect 
to the object. 
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Therefore 


du = Cd(oA) .< 3 > 


P 


T 


According to momentum equation, the net force across the vave front 
is equal to the change in momentum across the vave front 

• l 

• CA - d( + S'A = -(C-du)m 


or, 


Equating du yields. 


mdu 

du 


-d( A e ) 

- A 0 d<T ( 4 ) 

c 7>7 


From ( 2 ) 


d(yS>A) 

C2 


r 

- A. d<T 

o 

d(^5 A) 

*M£_2 
y ( 1-+6 > 

( 1 + f ) 2 d6~ 
y?o dT 


In the felastic region, the slope of the stress-strain curve is equal 
to the modulus of elasticity, giving 

C 2 = ( 1 + t ) 2 e 

which implies that the velocity of disturbance is depenent on strain. 
In Eulerian co-ordinate, the distance covered by the disturbance in 
time dt is, 

dx = ( u + C ) dt 

From the definition of total change in x 


dx 


= fb x ) dt + /3x \ da 
(dt/a (7a yt 


Equating the two equations yields: 


+ C = / dx \ da 

(writ 
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Therefore, the vave speed in terms of the Dagrangian co-ordinate system 



This vave speed is sometimes called the shift rate, =-g, in the solutions 
by the method of characteristics. . 

Therefore, 


.2 - 


“ m 2 


c 2 ' , 
(TTFY 


E_ in the elastic region, 




Po 

or 

g 

= ±fE 
Jp> 

Define 

the impact function <p 

as: 


b<p 

= 

or 

b«p 

* rb<f 



p a g 


From Newton’s lav, the momentum equation of an infinitesimal slice of an 
object at constant time can be derived as; 


or 


therefore 


+ A da - 5A 

° °3~a ° 

du - 1 

bt yO£> aT 

u - g 3q> 

"bt 2TT 



Aq da ) 3u 
C> t 



( 5 ) 


Differentiating (1) with respect to t yields, 

b fo X \ = £ u = ^ = l3q> 

<Tt^ a J *b a *9 t gbt 


therefore gbu - = 0 ••••••••••••••••••••••••( 6 ) 

b a ^ t 



Combining ( 5 ) and ( 6 ) yields, 

B (u-cp) + g5(u + cp) • =» 0 

St - "3a 

and 

3 (u + cp)-g2(u + cp) = 0 

d t "5a. 

Therefore ( u + <p ) r gS ( uv± <p ) = 0 

dt era 

But d(uj9') » <L ( u r <p ) dt + ^ ( u - cp ) da 

• dt .^a + 

d ( u - <p ) *=3 (u-®).^3 ( u T cp ) da 

dt + c>t ^5”a dt 


also 


± g 


da 

dt 


Therefore d(urep) = B ( u - cp ) + gS ( u $ cp ) » 

dt + ^a 

• IZT yhieh is the Characteristia equation. 

The equation establishes that, along the lines of slope + g 


( u + «p ) equals to a constant. 



B* Boundary and Initial Conditions. 


At the interface of any two materials, the forces and the 
particle velocities are equal: 



Material B 


And 


Material A 


In the elastic region: 

P « $A 

\ 

Therefore 


6r 


*kfk ea A A 

*A 


«B 

f B 


V/> £ A 

t PB/° B SB a b 

i 

» 

<pb « S 3 A u 
y^A £A Aa 


( 7 


< 8 


IP the structures are of the same material the relationship above 


reduces to 
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In this thesis the body of the specimen is initially at rest 
and undisturbed, and thus the field variables vanish identically 
ahead of the first vave front. 

If, at time t - 0,a continuous particle velocity (or stress) 

is applied at the boundary, discontinuity in cp and u will 
propagated along with the wave front. The characteristic equation 
can be conveniently used to compute the magnitudes of the propagating 
discontinuities. 



CHAPTER II 

• EXPERIMENTAL INVESTIGATIONS 

A. Experimental Arrangements 

Specimens are impacted by a projectile issuing from a 0.?5 inch 
inside -diameter air gun. The air gun consists of two main parts; a high 
pressure chamber and a barrel, 120 inches in length. It is mounted on 
a steel I- Beam and maintained at a horizontal level. Along the air 
gun barrel, vented holes are drilled to let out the air pressure behind 

. # i *\ ' ' * * 

the projectile to provide a constant impact velocity. The pressure 
chamber is pressurized by air from the air main. 

In the experiments, the compressed air can be released by 

opeming a solenoid valve operated by a 12 Volts D.C. battery. The 
■ *■> ■ -«. ‘ 
projectiles made from machined aluminum rods of various lengths are 

loaded through the open end of the gun barrel. Two adjustable teflon 

supports are mounted on the I- Beam close to the open end of the air 

gun. Teflon was used to reduce the friction between the support and the 

specimen after it was hit by the projectile so that free impact could 

be achieved. ( Fig. 1 ) 

The specimens were machined from 10 inches long, 0. 75 inch diameter 
aluminum rods to the proper configurations ( Figs. 2 and 3 )• Four 0.25 
inch long dynamic strain gages were mounted on each specimen as sensing 
devices. The gages were mounted in. diametrically opposite pairs to cancel 
any stress due to bending effects. 


n 



Schematic Diagram of Experimental. Set-up 
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Grooved Specimen and the Experimental Parameters 



B* Ehqperimental Procedure 

Each pair of strain gages was connected to the opposite 
arms of the Wheatstone circuit of a bridge amplifier together with 
resistors as dummy gages. For discussion purposes, the location 
of the first pair of strain -gages from the impact end of the specimen 
is called station 1, that of the downstream pair is called 
station 2, respectively. The bridge amplifiers were calibrated at 
the same gain for the consistancy of results. Every specimen was 
grounded to eliminate noise which might be introduced through 
the gages. 

The output signals from the bridge amplifiers were 
connected to a dual beam oscilloscope. The signals from stations 1 
and 2 were given by the upper beam and the lower beam^respectively. 
The oscilloscope was set at internal trigger and single sweep mode 
using the signal from station 1 as the triggering signal. 

Soon after the projectile hit the specimen, the first 
stress pulse which passes through gages at station 1 is transmitted 
through the bridge amplifier to the oscilloscope. The rising part 
of this pulse triggers the single sweep .setting and gives the shape 
Of the pulse. At the same time the lower beam which represents the 
signal from station 2 traces along the reference line until the 
gages "feel” the first stress pulse passing through. 
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C. Experimental Programs. 

To study the configuration effects in stress wave propagation, 
specimens were made with strain gages mounted at various positions as 
shown in Figs. 2 and 3» The tolerance of machining of the specimens was 
+ 0.005". In Ike grooved specimen (Fig. 2), the length of the groove in 
one set of experiments was changed from 0.125 H to 2.937"with the diameter 
of the groove held at 0«531"» In another set of experiments, the diameter 
of the groove was changed from 0.40" to 0.70" with the length of the groove 
held at 2.00". 

To study the effects of threads in stress wave propagation, 
experiments were conducted using specimens with plain holes, threaded 
holes and threaded holes with screws. The parameters were the depth and 
the diameter. The threads were American Standard and the sizes corresponding 
to: the sizes of the holes being used in the investigations. 

In: all above experiments, three different lengths of projectile were 
also used. i.ir:.: : -r . 

'To study the stress decay in regions close to a section of 
discontinuity, specimens similar to those in Fig. 2 were made. Three sets of 
experiments were conducted, in the first two sets the diameters of the 
grooves were held at 0.531" and 0.45" respectively while the length of the 
groove was varied from 0.125" to 2.50". In the third set of experiment, 
the distance from station 2 to the downstream edge of the groove was 
held at 0.10" while the diameter was changed from 0.40" to 0.70". 

Experimental programs with the grooved specimen and the hole 
specimen are summarized in Tables 1 and 2. 


I 
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TABLE 1 

Experimental Program with Grooved Specimen 


Length of Groove^ A ( in. ) 

Diameter of Groove-^ (in.) 

Projectile length (in< 

0.125 

0. 760 

1.00 

0.250 

0. 600 

1.50 

0.500 

0.575 

2.00 

1.000 

0.550 


1.500 

0.531 


2.000 

0.525 


2.500 

0.500 


2.750 

0.475 


" 2 . 8?5 

0.450 


2.937 

0.400 



D© = 0.531", C = 2.000" A = 2.00", C =2.00" 

D o * 0 * 4 50", C = 3.000" 


A 


2.00", C = 3.00" 
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CHAPTER HI 


RESULTS, DISCUSSIONS AND CONCLUSIONS 


A# Comparison of Experimental Results to Numerical Solutions 

A computer written by Tang ( Ref. Appendix I) was used to obtain - 
results to numerically check against those of the experiments. The 10" 

.long .specimen was divided into 100 imaginary layers of 0.1" in thickness. 
Because the computer program is based on the calculation along grid networks 
the smaller the characteristic grid, the more accurate the results willrbe. 
For an input pulse of duration 20//. sec and a wave speed of 0.2 inch^-sec 
in aluminum , the input pulse will be sampled into 20 sample periods of 
l^sec. each. This sample period yields acceptable results. 

At a boundary of discontinuity, (Fig. 5), the stress magnitudes of 
the reflected and the transmitted waves are related to the magnitude of 
the incident wave by the formulas. 


<Tr 

Vhere ^ 

<Tt 

. So 

A o 

A t 


= C( At - A c ) 



= Reflected Stress 
= Transmitted Stress 
= Incident Stress 

= Original Area 


Transmitted Area 
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.Fig* 5 Stress VJave at a Section of Discontinuity 

The. stress and particle velocity at the free end must 
. .satisfy ‘the end conditions; In this case, the free surface — 
can not be subjected to any longitudinal stress, therefore the 
©tresses of the incident and reflected waves must add up algebraically 
to aero. The sense of the stress in the reflected wave is the •*, 
opposite of the incident wave and the magnitudes are the same. 

The particle velocities for the reflected waves add up to double 
the particle velocity in the incident wave. (Fig 6) 



• Fig. 6 Stress Wave at a Free End 



Numerical computations were done for stress wave propagation in 
grooved specimen and in a. -hole .specimen. In the computations, the input 
stress pulse given by the gages at station 1. was approximated by a sine 
wave. The comparison -of the numerical computations to the experimental 
results is very good as can be seen from Figs. 7 and 8. 



Numerical 



Fig* 7 Comparison of Experimental and Numerical Results in Grooved Specimen, 
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B. Configuration Effects. 

1* Effects of Grooves, Holes, Threads, Threaded Holes and Threaded Holes 
with Screws. 

As the stress wave propagates through the discontinuity, part 

of the wave is reflected and part transmitted; their magnitudes are 

dependent upon the area ratio at the boundary. This phenomenon repeats 

itself at the next discontinuity but with a time difference. The 

stress at any point in the specimen at any time is the algebraic sum 

of the magnitudes of the waves arriving at that point at that time. 

A hole in a specimen constitutes a reduction in area in that 

section of the specimen. For one dimensional analysis and calculation 

of stress waves propagation, the hole section can be treated as a solid 

rod with a cross- sectional area A„, smaller than that of the specimen. 

r 

A q = Cross-sectional Area 
of Specimen 

, ^A^ = Cross-sectional Area 

of Hole 

A r = Cross-sectional Area 
of Solid Rod 
L = Depth of Hole 

© 

Fig. 9 Modification of Hole Parameters in Calculation. 
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For the threaded holes, American Standard taps of diameters equal to 
tuthose of the plain holes were used. 

In this case, the cross-sectional area is smaller than that of the 
plain hole due to the threads. Since the depth of the threads- is small 
compared to the diameter ( see Table 3 )» the overall effects to stress 
wave propagation are expected to be small. In calculations, an effective 
area of ( A Q - A^ + 1 A^ ) where A^ is the annular area of the threads i 
can conveniently be used. 

TABLE 3 

American Standard Threads 


Thread/inch 

Outside Diameter (inch) 

Thread Depth (inch) 

20 

0.25 

0.019 

16 

0.375 

0.037 

' 13 

0.50 

0.040 


\3hen a screw is tighten in the threaded hole, the specimen 
theoretically will act as a solid rod. However, due to machining 
errors there are gaps between the threads and the screw which can 
be considered as a small hole in the -rod. The magnitude of the 
stress is expected to be higher than that of the solid rod. 

Experimental results of stress wave propagation in rods with hole 
threaded hole and threaded hole with screw are presented in Fig. 10. 
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The comparison of the maximum stresses rj nax at station 2 to the 
input stress pulse (S' , for the cases of plain hole, threaded hole and 

threaded hole with screw are plotted in Figs. 11 and 12. 

From the plots, it can be seen that as the diameter of the hole 
becomes larger, giving smaller annular area for stress waves propagation, 
the stress peaks gradually increase. For the case of threaded hole with 
screw, the peak values approach a certain constant value slightly larger 


than that of the solid rod 
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2. Effects of the Length of Projectile 

fhe stress wave is generated by collision of the projectile to 
the specimen. By the method of characteristics, it can be easily shovn 
that the duration of the pulse is equal to twice the time elapsed T in 
stress propagating through the projectile’s length. The field variables 
tt and cp for the projectile and the specimen in Fig 13 show that contact 
will be maintained for a duration of 2T.vhen = V . After this 

duration, there will be no contact since u^ra.u^.- 



From the experiments, the duration of the applied pulse is 
slightly longer than theoretical result because the alignment in 
colliding is not perfect. 



3 * 


The experimental results of the pulse durations for the 
-projectiles of length 1", 1.5" and 2" are shown in Fig 14. 


a) 1" projectile pulse duration 


"b) 1.5" projectile pulse duration 



Fig. 14 Experimental Results Showing Change of Pulse -Duration 
Due to Various Projectile Length 
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C. The Stress Decay in Regions Close to a Section of Discontinuity 

The behavior of stress wave propagation in discontinuous structures 
resembles the flow of a fluid through pipes of unequal diameters (Fig. 15). 
It will be some distance downstreams from the discontinuous boundary before 
the stress resumes its predicted proper magnitude. 



Fig. 15 Stress Flow Analogy at A Section of Discontinuity 

To investigate the phenomenon of stress decay in the vicinity of 
an area of discontinuity, grooved specimens with step Y 0 (Fig. 15) were made 
with strain gages mounted at a distance, Xq, downstream from the section 
of discontinuity. The stress decay can be -detected by comparing the first 
peak stress at station 2 to that obtained from the theoretical calculations. 
This first peak stress is the stress felt by the strain gages at station 2 
as the input stress wave first swept by. For a sufficiently large value 
of Xo , the experimental and the theoretical peak stresses should be equal. 
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The program of investigation has been outlined in .the section 

of Experimental Programs. Briefly, vith Y =0.109", X was varied from 

o o 

0.15" to 1*00". The result of the comparison, of the peak stress, 

to the input stress (f 0 , was plotted in Fig. 16. It can be seen that 

a distance X = Y from the section of discontinuity is reauired for 
o o 

the stress to completely resume its predicted magnitude. 

Additional experiments were conducted with X Q = 0*10" while Y 0 

was varied from 0.^0" to 0.70". The ratio of to (T 0 was plotted in 

Fig. 17. The result indicates that the gages at station 2 detect very 

little stress for value of Y c ^1.5 

X o 

A sample result of the experiment ( X a - 0.125, Y^= 0.109" Yq^D.87) 

*o 

as compared to the numerical solution is presented in Fig. 18. 






ssaaq-g 


i Close to a Section of Discontinuity 
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D. Conclusions and Recommendations 

From the analysis of the experimental results, it is concluded 

that • 

1) It is possible to apply one dimensional wave theory to obtain 
stress information for discontinuous structures. ( A computer program 
.. . has been written for the purpose of obtaining such information.) 

Z) Small discontinuity in structures such as threaded hole, 
threaded hole with screw has little effect in stress wave propagation. 
In the analysis plain holes using an average area can conveniently 
be used in place of the threaded holes. For design purposes, structures 
consisting of threaded holes and threaded holes with screws can be to 
conservatively approximated using dimensions of the plain holes. 

•• 3) As in the cases of fluid flow through the pipe sections, 
streamlining will not be reached until some distance downstream. 

Stress wave propagation shows the same behavior. Experimental results 
indicate .that to obtain the proper stress magnitude, strain gages 
should be placed at a distance from the section of discontinuity at 
least equal to times the length of the 'step. 
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APPENDIX I 


NUMERICAL PROCEDURE AND COMPUTER PROGRAM ( Ref. 4 ) 

Consider the grid system of the characteristic lines in the 
a-t plane as shown in Fig 19« The ordinate a= a^ represents 
the boundary of the specimen, where the input s(external disturbances) 
are applied. The numerical procedure involving stepwise integration 
along the characteristics of u± cp = constant along the line of 
..slope T g is employed to compute the response for various inputs. 

o' 

The stresses and the particle velocity are calculated at all 
grid points. 



Fig 19 Characteristic Network 



Since the external input has a finite rise time, p and u 
are identically zero along the wave front A B G the grid points 
are devided into 3 groups 

- interior points 

- boundary points 

- interface points 

Boundary Points 

-Along the boundaries at a = a^ and a - a Q one of 
the two variables p and u is perscribed as a function of time. 
The remaining variable may be determined by representing the 
characteristic equation in finite difference form along the 
characteristic curve referring to Fig 19» since the field variable 
vanished at A and B for the interior boundary a = a^ , only 
the equation along g - (C to 3 is used). For the exterior 
boundary a = a Q » only the equation along g + (E to F) is 
employed. 

3 Interior Points 

The field variables vanish at B and G, and they have been 
computed at the boundary point C, the value of P and u at the 
grid point H can then be computed by employing the characteristic 
equations in finite difference form along CH (g+ ) and GH (g-) 


40 


The computation at other interior points proceeds step by step 
procedure. Thus, to compute the field variables at S , we employ 
the values at P and R, which have previously been computed, and 
we-use the characteristic finite difference equations along g- 
and g +, respectively. 



the points and » "which are located on different sides of 
an interface, but very close to the interface, as shown in Fig 20. 
There are altogether four unknown field variables at the points 
Nj—and • Above the interface, the characteristic equation 
.along can be utilized. Underneath the interface, we can use 

the characteristic equation along LN^. Two additional equations, 
adding to a total of four equations, are obtained from the interface 
conditions (7) and (g), which relate the field variables at 
and 



The four field variables at the interface can thus be computed 
by solving four algebraic equations simultaneously. 

The computations described above were carried out on a 
digital computer. It is apparent that the smaller the characteristic 
grid, the more accurate the numerical results will be. 

Care must be exercised to ensure that the g + characteristic 
lines reaching an interface meet up with the g - characteristics 
of the next layer which carry the reflected waves from the next 
interface. This can be accomplished by relating the increments 
/i a in the ith and (i+l)th layer by the relation 

' e i ^ H+i* a i+ i 

vhere and g^ + ^ are the dilatational wave speed In the ith 
and (i+l)th layers. 
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•RUN T0D5AR, 306-01-122, YANG' 

•FOR* IS MAIN 

C ( N ) = C(I)/C(I + 1J RATIO OF WAVE VELOCITY 
E(N)= EII1/S0 RATIO OF MODULING TO INPUT STRESS 

K IS THE NO OF LAYERS . 

KK(N) is THP LAST GRID PT. OF EACH LAYER 
V ( N ) IS TH<^ POISSON RATIO 
S ( 2 * I ) IS THE RADIAL STRESS 
U ( 2 , I ) IS THE PARTICLE VELOCITY 

JT.JE IDENTIFIES THE TRIANGULAR INPUT 
L IS THE TIME THAT PROGRAM ENDS 
DIMENSION U( 2 , 1000 ) , s < 2 ,1000 ) 

DIMENSION KK ( 3 ) , H ( 2 ) 

REAL J1,JT 
REAP 1 ,L,K* JT 

1 FORMAT ( 2X,2I4,F5. 1) 

READ 2 » H ( 1.) *H( 2 ) 

2 FORM \T ( 2X » 2F5 « 3 ) 

REAt 3 » KK ( 1 ) , KK ( 2 ) » KK < 3 ) 

3 FORf* \T ( 2X * 3 I 4 ) 

Kl = * < { K ) 

DO 60 1=1* K1 
S(1»I)=0. 

U(1 ,1 )=0. ' • 

6 0-C0NT4NUE 5 

DO 500 J = 2 » L 
J1=J-1 

IF ( Jl-JT) 75,75,78 
75 R=J1/JT 

5(2,1 )=SIN(3.1417*R) 

GO TO 90 
78 S( 2 , 1 ) =0. 

90 U ( 2 , IT -U ( 1 » 2 ) +S ( 1 » 2 ) -S ( 2 ♦ 1 ) i 

DO 400 N= 1 , K t ! 

M=KK(N)-1 

IF (N-l) 124,124,125 - 

124 NI =2 ' 

GO TO 126 

125 NI=2+KK(N-1) 

— 1 2 6 DO 130 I = N I , M ' 

S(2,I) = (U(1,I + 1)+S(1,I + 1 ) — U C 2 , I -l)+S(2,I-l))/2. 

U(2,I ) = ( U ( 1 , I +1 ) +S ( 1 , I + 1 ) +U ( 2 , I -1 ) -S ( 2 , I— 1 ) ) /2 • 

130 CONTINUE 

IF (N-K) 140,420,420 
140 KMMM=KK(N) 



MM=KMMM+1 
MMM = K K ( N ) + 2 ' 

S(2 sXMMM) = (U ( 1 »MMM)+s( 1 2 »M)+S( 2 »M > ) 

S ( 2 ! CMMM } = S ( 2 , K.MMM ) / ( 1 • +H (Ml) 

U 12 * CMM M) = (U( l,MMM)+c( 1 ,MMm ) +H ( N ) * ( U ( 2 tM ) -S ( 2 » M ) ) J/ll.+HCN) ) 
S(2? 4 M)=h(N)*S(2» KMMM ) 

U(2,MM)=!J(2,<mMM) 

400 CONTINUE 

420 KV„MN = K:<( K) • 

. S(2»KMMN)=0. 

U ( 2 » KMMN ) = U ( 2*M)— 5(2»M) 

250 PRINT 255,J>S(2,1),S{2,21),S(2»82)>S(2»83) 

255 F0RM4T (2X» I4.4F10.6) 

DO 490 1=1, KMMN 
S(1 * I ) =S ( 2 » I ) 

U(1 »I)=U(2»T ) 

490 CONTINUE • 

500 CONTINUE . . 

STOP 

END 
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LIST OF INSTRUMENTATIONS 


ITEM 


. SPECIFICATION 


Strain Gage 


Micro-Measurement, Type ED-DY-125AD-350 
Resistance 3 50.0 0.3 $ ohms* 

Gage Factor 3-2 1 2.C^ at 75° F 
Gage_Length 0.25 inch : — — 


Bridge Amplifier 


Ellis-Yishey Bridge Amplifier Model 
BAM-1B, Cutoff frequency 130,000 C/s 


Oscilloscope 


Textronix 'type 551 Dual Beam Oscilloscope 
vith 1A5 and CA Plug-in units 

Alpha Shielded ’.■fire 


\Iire 
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ABSTRACT 


Title of Thesis: Radial Elastic-Plastic Stress Wave Propagation 

with an Unloading Shock in a Circular Disk by 
the Method of Characteristics and Successive 
Elastic Approximations 

Donald W. Seaton, Doctor of Philosophy, 1972 

Thesis directed by: Dr. Jackson C. S. Yang 

Dr. Arnold E. Seigel 

The purpose of this work was to investigate the problem of 

radial elastic-plastic stress wave propagation in a thin, finite, work 

hardening disk using the method of characteristics and the method of 

successive elastic approximations. In this method the characteristic 

equations, derived for the elastic case, are applied to the plastic 

regime by allowing the modulus of elasticity E and the wave speed C 

% 

appearing in the equations to assume values corresponding to the tangent 

t 

modulus of the plastic stress-strain curve. The value of Poisson's 
Ratio, v, is taken as 1/2 in the plastic range. 

Three characteristic equations were derived from the governing 
conservation and elastic constituent equations as outlined in a paper 
by Chou and Koenig. They were used subsequently, along with the Prandtl- 
-Reuss incremental plasticity theory and the method of successive elastic 
solutions , to calculate the problem variables at the nodes of the 
characteristic network which was generated in a step-by-step fashion 
as the solution progressed. 

A modest experimental program consisting of two experimental 
techniques was conducted to check the accuracy of the proposed analytical 
method. A circular disk was used with a hole at the center to accommodate 


a pressure pulse applied either hydraulically or explosively to the 
Inner surface. The hydraulic method produced poor results and was 
discarded in favor of the explosive technique. 

The disk material was 1100-0 aluminum which exhibits a 0.2 per- 
cent offset yield strength of 5000 psi and a proportional limit of 2500 
psi. Strain-rate effects were not included in the theory. 

Two 350 ohm Micro -Measurement strain gages were used in addi- 
tion to a TEKTRONIX Model 551 dual beam oscilloscope to obtain radial 
and circumferential strain time histories at a radial distance of 0.75 
inches on a 4 inch O.D. specimen. 

" * — A" Piezotronics Model 109A ballistics pressure transducer was 
used with one channel of a second Model 551 TEKTRONIX dual beam oscillo- 
scope to record the applied pressure pulse. A computer program in- 
corporating the analytical procedures described previously was used to 

check the strain histories measured on the specimens. In addition, the 

% 

theoretical solutions for a ramp followed by a steady state pressure were 
run as approximate checks against the known, closed form, static elastic 
and elastic-plastic solutions for a thick-walled cylinder. 

Reasonable agreement between the theoretical and experimental 
results supports the hypothesis that this method can be used to solve 
problems of wave propagation in systems subjected to combined states of 


stress 
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SECTION I 

' > 

INTRODUCTION 

The theoretical study of plastic stress wave propagation in 
solids first received serious attention in 1941 with the pioneering 
work of Von Karman : and Duwez, later published in 1950 [13, and the 
concurrent work by . G. I. Taylor in England [2] and Rakhmatulin in the 
-Soviet Union [3] . . 

Out of this work came a strain rate independent theory for the 
propagation of one-dimensional, plastic waves. However, as mentioned 
in [1], strain rate effects were observed in experimental results for 
which the rate of ^.training was reasonably above that encountered in 
Static or quasi-static loading. Accordingly, later Investigators such 
as Malvern [43 modified the basic governing equations to include the 
effects of strain irate which resulted in improved agreement between 
theoretical predictions and experimental results. Active research has 
continued in an effort to develop more accurate experimental techniques 
to test the validity of various theoretical formulations as evidenced, 
for example, by the work of Bell [53 • 

Considerable work is also being done on the propagation of com- 
bined states of stress in one or more spatial dimensions such as that 

% 

presented by Clifton [63, Ziv [73, Chou and Koenig [ 83 , Fyfe and Swift 
[93, and Lawrence [103 • 

It is the purpose of the present work to develop a method of 
solving the problem of radial elastic-plastic stress wave propagation 
with an unloading shock in a thin, finite, work hardening circular 
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disk, by combining the method of characteristics and successive elastic 
approximations utilizing the Prandtl-Reuss incremental plasticity 
theory. A computer program utilizing this method and included as 
Appendix D, calculates the time history of the problem variables. 

By using final values of the appropriate variables as initial values 
for a second run, the program can also simulate the case of repetitive 
loadings on a single specimen. 

The results of a modest experimental program confirming the 
theoretical solution is also presented. 


SECTION II 


PROBLEM DESCRIPTION 

As mentioned in Section I , the present work deals with the 
theoretical and experimental analysis of the propagation of elastic- 
plastic stress waves in a finite, circular disk with a central hole in 
which is initiated a symmetric, radial pressure pulse. The outer 
boundary of the disk is considered free. 

The characteristic equations governing the propagation of 

elastic stresses and strains are derived in detail in Section III along 

with the relationships describing the method of successive elastic 

approximations. A thorough discussion of the computational procedures 

is given in Section IV. It is sufficient to say here that as the 

material is loaded into the plastic range, the elastic characteristic 

% 

equations are applied in an incremental fashion to each of the character 

I 

istic network nodes which are generated one-by-one as the solution pro- 
gresses. During this process the modulus of elasticity E, Poisson's 
Ratio, Vj and wave speed C, appearing in the elastic equations are ad- 
justed to values corresponding to the local tangent modulus of the 
plastic stress-strain curve. The value of Poisson's Ratio, v , is taken 
as 1/2 in the plastic range. 

In addition to the propagation of elastic-plastic character- 
istics, the problem of an unloading shock wave is also treated. This 
occurs whenever the equivalent plastic stress at any radial coordinate 
decreases with time. The equations describing the unloading wave 
phenomena are also derived in Section III. It is assumed during 
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unloading that reverse yielding does not occur. 

During experimental testing two methods have been used to 
initiate the pressure pulse. The first consists of applying the pres- 
sure hydraulically by dropping a weight on a piston which compresses a 
quantity of water constrained in the central specimen hole. This pro- 
cedure, however, yields unsatisfactory results. The second method 
uses a fuse ignited explosive charge and has been found to yield more 
realistic data. 

For the sake of completeness the experimental fixture used in 
the hydraulic method and typical test measurements are included although 
computer solutions attempting to simulate these measurements are not 
presented, but instead are restricted to the- results simulating the ex- 
plosive tests. 

As a result of the applied pulse, a stress wave propagates 

radially outward to the free surface where it reflects and propagates 

inward to repeat the cycle. The amplitude of the pulse has been chosesn 

» 

so that the material behind the wave fronts is eventually stressed into 
the plastic range throughout a central, annular region of the disk. The 
state of stress is one of plane stress wherein only the radial and cir- 
cumferential components are significant. There are, however, radial, 
circumferential, and axial components of strain. 

The strain history is measured at one radial location with two 
"350 ohm Micro-Measurement strain gages mounted 90® from one another on 
the flat surface of the disk. The sensitive axis of one gage is oriented 
along the radial direction while the axis of the remaining one is posi- 
tioned along the tangential or circumferential direction. (See Figure 
lb.) The gage outputs are recorded on a Type 551 TEKTRONIX dual beam 



Oscilloscope and polaroid film. 

-The input pressure pulse signature is measured with a 
-Piezotronics quartz transducer, Model 109A, mounted in the support 
plate with the input surface centered immediately below the central 
"tustd in the specimen. (See Figure la.) The output of the trans- 
-ducer is also recorded on one channel of a Type 551 TEKTRONIX dual 
beam oscilloscope and on polaroid film. 

A more complete description of the experimental procedure and a 
detailed sketch of the apparatus can be found in Section V and Figures 
la - le respectively. 
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e Z1 - • ■ - 

. SECTION III 

DEVELOPMENT OF GOVERNING EQUATIONS 


- - - If we consider an infinitesimal element removed from a sector 

of a disk and subjected to a circumferentially symmetrical state of 
plane stress as shown in Figure (2), we may easily show that the equi- 
librium equation in the radial direction is 



dv 

at 


3-1 


Here, and in what follows, we are referring to polar, Lagrangian type 
coordinates where the origin of the coordinate system Is fixed at the 
center of the disk. There is no equilibrium equation in the 0- 
direction since we are assuming circumferential symmetry. Furthermore, 


since plane stress is being considered there are no stress components 
in the Z- direction. Two other equations which are at our disposal 
include the elastic constitutive equations; i.e., the stress equations 
in the r, 0, and z directions. It will be seen, however, that only two 
of these are required to provide three equations in three unknowns; 

* a , a B , and the particle velocity v. 

The constitutive equations are 
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3-4 

time 

3-5 
3-6 

Equations 3-1, 3-5, and 3-6 are sufficient to derive the 
characteristic relationships which consist of a set of differential 
equations relating the dependent variables a , a Q , and v to the inde- 
pendent. variables r and t along the characteristic lines in the r-t 
plane. These are lines along which weak disturbances propagate and 
across which the dependent variables' are continuous but the derivatives 
of the variables may be discontinuous. 

There is more than one way to derive these equations, (see [6} 
for instance) but in view of the above definition, it is perhaps easiest 
« to solve 3-1, 3-5, and 3-6 directly for the derivatives by Cramer's rule 
and then to impose conditions which will cause the derivatives to become 
undefined. 

Accordingly, we rearrange the equations in the following form 
and also include the general relationships for the differential changes 

-Of o , a., and v as a function of infinitesimal changes in r and t. 

r . y 


** * Hr ■ r [- ” <a r + v! 

o 

Differentiating the first two of these with respect to 

&u 

and noting that the particle velocity v = — gives. 


3v 

1 

r 3a r 


3r 

~ E 

Ut 

” St J 

v 

_ 1 

r aa e 

3a -i 
„ — 1 

r 

= E 

L3t 

v at J 


Characteristic Equations 



bo. 


r 

bx 


bv 
* p b£ 


bo. 

Tr 


, bo 

r 

E bt 
bo 

. H __L 
E St 


bo 

1 


v ^ CT 6 _ Sv 

E St " Sr 


+ I^i 

E St 


g fl l q r 
r 

0 


v 

r 


dt 


So ft Sa fi 

dr + -r- dt 
Sr St 


Sv 


Sv 


* dr + rr* dt = 
Sr St 


bo- 


do 


da. 


dv 
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We may solve the above system for — ~ , for example, by Cramer's 

or 

rule. The result is a fraction consisting of a determinant in both 
numerator and denominator. If the denominator is set equal to zero and 
expanded, we get a relationship describing the slopes, or directions, 
of the characteristics in the r-t plane. Similarly, if we set the 
determinant in the numerator equal to zero we obtain the differential 
equations relating the unknown variables along the characteristics. 

(See Appendix B.) The same relationship would be obtained if we chose 
to solve for any of the other derivatives. - 

The results of these computations show that the characteristic 
directions are given by the following three equations 




3-8 


dr 

dt 


- c 


3-9 
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dr « 0 


3-10 


Along the first of t hese direct ions . 


dcr r * pc dv = £- (a r - a Q ) + pvcvj 


3-11 


Along the second. 


d<J r + pc dv 


■[ 


' (a r ’ V 


- pvcv 


']* 


3-12 


and finally, along the third. 


ia T ’ » [% ' E r dt ] 


3-13 


In these equations E, c, and u are to be regarded as variable 
quantities the values of which depend upon whether the state of stress 
ritt .the material is in the elastic or plastic regime. 


Plasticity Relations 

% 

Equations 3-1 to 3-13 are sufficient to define the state of 
stress and strain in a perfectly elastic, homogeneous, isotropic. 

Strain rate independent material at constant temperature. However, 
vhen the stresses become sufficiently high a yield point is reached 
and the material behaves plastically thereafter or until unloading 

« 

occurs. In the latter case unloading proceeds along a curve parallel 
- -to the original elastic curve resulting in a permanent set when the 
stresses reach zero. Upon reloading, the material remains elastic 
until the stresses again become equal to the maximum value previously 
attained. (See Figure 3.) It is also assumed in what follows that the 
stress-strain curve in compression is identical to that in tension. 
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A detailed development and discussion of the various plasticity 
theories in common use is given in[H3as well as in a number of other 
excellent references. Therefore, only a brief outline will be given 
of the theory which is used here. 

In the case of uniaxial loading, a material specimen will 
behave elastically only until the applied stress reaches the yield 
point £Tq. Beyond this value plastic deformation occurs and the 
relationship between the stress and strain is no longer linear; i.e., 
yielding is said to occur when a applied 2: Oq. When a multi-axial 
stress state prevails one of several possible yield criteria is used 
to define the onset of yielding. In this work the Von Mises or Dis- 
tortion Energy criterion for biaxial stress is used in which yielding 
is assumed to occur when 


r 2 2 

CT + o* - a Ca 2: C„ 
L r 0 r 0J 0 


3-14 


This function, however, is also the definition of the equivalent 

stress, a , which represents the ’’effective stress” on the material, 
eq 

Thus, yielding occurs when 


a *= 
eq 0 

When < 7 ^ becomes greater than o^, the functional relationship 
between the stresses and strains is non-linear. The extension of 
3-2, 3-3, and 3-4 to include the plastic strains can.be written 


e„ *= “■ {~cr - vcr a 

1 + e P + de P 

3-15 

r E L r 9. 

Jr r 

’ c e ” e [°e " ” CT r. 

] + Cq P + dc Q P 

3-16 
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«. ■ E L" v ( °r + ?e } ] 


8'J + e l + d£ z 


3-17 


These equations suggest that the strains are computed from the 

applied loads in an incremental fashion. This is done with the aid of 
the Prandtl-Reuss equations which are given in general differential form 

as ' * ; 


de ij “ S ij dX 


3-18 


where is the deviatoric stress tensor and represents the stress com- 
ponents which are in excess of the mean hydrostatic stress which. pro- 
duces no yielding. For the case at hand, this can be written 


ij 


a - cr 
r m 


0 

0 


Ov - a 
6 m 


■a - a 
z m 


Since the mean stress is defined as 


m 


5 [° r + °e + ° J > 


*the deviatoric stress tensor finally becomes 


ij 


2o r - a Q 


3 

0 


2a, 


a 

. z 


* a 


3 

0 


2cT z " a r ‘ a e 


3-19 
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By invoking the yield criterion it is shown in [HI that the constant 


d\ *= 3/2 


where dtp is the equivalent plastic strain increment and is given by 


'P 3 


[<de* - de^) 2 + (de^ - de^) 2 + (de* - de*) 2 ] 3-20 


Substituting this relationship and 3-19 into 3-18 gives the following 
expressions for the plastic strain increments in finite difference form 


with a — 0 . 
z 


. p Ac p r\ 1 
Le r c iTT L 2a r " CT eJ 


3-21 


P Ac P 

Ae e * 2a a 
eq 


[ 2CT e - a rJ 


3-22 


4 P , P. . P 
Ae = - Ae - Ae Q 
z r o 


3-23 


Equation 3-23 follows from the fact that the first order approxi* 
mat ion for the relative change in volume. of - a strained element is. 


= e + e Q + c 
v„ r e 2 


which in the plastic range is neglected in comparison with the deforma- 
tions. Furthermore, it may also be shown that the relative volume 
change is given by 
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which implies that v - 1/2 when the change in volume can be neglected, 
i.e., in the plastic regime. 

Equations 3-21, 3-22, and 3-23 are the so-called Prandtl-Reuss 
flow rules for the incremental plasticity theory. 

- - -Only one additional set of equations is required to define the 
stress-strain time histories. These are the unloading shock relation- 
ships to be discussed next. 


* 


Unloading Shock 

If at any radial distance r the material is stressed beyond the 
yield point so that ^ SIGMAX, the maximum equivalent stress pre- 
viously attained, and subsequently decreases with time, then unloading 
is said to occur as depicted in Figure (3). Since unloading proceeds 
along an elastic curve, such a "disturbance" propagates with the elastic 
wave speed c q . It therefore overtakes the slower plastic waves and may 
prevent them from reaching the outer boundary or may itself be prevented 
from reaching the boundary. The method of determining this behavior 
will be discussed later in this section. 


. As previously stated, the characteristic equations given by 3-8 
to 3-13 are applicable to problems in which the variables are continuous 
but the derivatives of the variables may be discontinuous. In the case 
of an unloading shock the variables themselves are discontinuous across 
the shock and are therefore not governed by equations 3-8 to 3-13. The 
shock equations, or equations describing the propagation of a discon- 
tinuity may be derived however, by employing equations 3-11, 3-12, and 
•3-13 in a limiting, process. The following closely parallels the 
derivation given in [8], originally credited to Leonard and Bud'iansky 
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v aiwj Jahsman. 

. Consider an infinitesimally thin shock wave, or discontinuity, 
propagating along a 4- c q direction in the r-t plane as shown in, 

Figure (4). 

Xf we. also consider, as shown in Figure (4), a - c q character- 
istic intersecting the shock discontinuity, we can write equation 3-12 
• between two points A and B on either side of the shock with the result 

that 

B 

(a iS - a rA> + l )c o (v B * V * J [- <a r ' V ‘ " ucv ] T 3-2 

A 

As B is allowed to approach A, or as dr •* 0,.the integrand contains 
bounded values of a , a Q , and v so that the right-hand side vanishes. 
Also, a _ - a . is finite but different from zero and represents the 

Tn TA 

Jump in at point A across the shock propagating along c . If we de- 

% 

fine 



with analogous expressions for the other variables, equation 3-23 can 
be written finally ... 

6a r + pc Q $v =0 3-25 

In a similar fashion, if we had considered the shock to be 
propagating along a -.c line and instead of using equation 3-12 we had 

c r ■“ 7 W 

written equation 3-11 along a + c^ characteristic intersecting the shock, 
equation 3-25 would be replaced by 

c-> . &° r “ pc &v = 0 - • 3-26 
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Finally, If equation 3-13 is applied between points A and D in 


Figure (4), we get 


U 

°rA v L CT 8D CT 6A J E r dt ] 


As dt -* 0, the integral approaches 0 and we are left with 


CT rD ” C rA 


■K«i 


0D ' CT 0A. 


and since in the limit cr _ cr _ cr . we have finally the static re- 

rD rB rA 

lation at a given value of r. • 


to r " 7 6a e 


3-27 


Equations 3-25, 3-26, and 3-27 relate the jumps in the variables 
a , Oa» and v. at a point along a discontinuity in the r-t plane. However, 

T V 

% 

we would also like to know how the amplitudes of these jumps vary as they 

* 

propagate outward or inward along + c q or - c q respectively. This is 
done for propagation along the + c q direction, for example, by writing 
- equation 3-11 along the right hand (upstream) edge of the shock while 
applying equations 3-25 and 3-27 across the shock. This provides three 
* equations in the three unknowns a , o 0 , and v at a desired point along 
"the upstream edge of the shock. The finite difference representation of 

• . • c. w . . . • 

these equations will be presented in Section IV. 


vr _ .. - . . Boundary Conditions 

• There are two boundary conditions in this problem corresponding 
to conditions on the inner and outer surfaces, or edges, of the disk. 
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The 'inner surface experiences an applied pressure as a function of 

time which, of course, is equivalent to radial stress as a function 

of time. The outer surface is stress free and hence a =0. 

r 
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SECTION IV 

NUMERICAL ANALYSIS PROCEDURES 

Equations 3-11, 3-12, and 3-13 are the differential equations 

which will be used to compute the quantities a , a. , and v in a stepwise 

r 0 

fashion along the characteristics in the continuous field region. 

Consider a typical r-t plane as depicted in Figure (5) with a 

plot of the input pressure (c ) shown directly below. We divide this 

plane into a number of lines of constant r which shall be designated by 

the index I running from I = 1 to I = NA. (In this analysis A stands 

for segment A. In the future, if one should desire to consider several 

concentric disks or rings these may be conveniently labeled B, C, D, etc.) 

The + c characteristics beginning at the 1=1 line and propagating 

(either elastically or plastically) to the I = NA line are denoted by the 

% 

index J and cover the range from 1 to problem termination or the limit JLIM. 

t 

By using equations 3-8 and 3-9 in addition to the three differ- 
ential, equations , we may construct the characteristic network of Figure 
(5) point -by-point in an incremental fashion. A more detailed description 
of this procedure will be given presently. 


Finite Difference Equations 

. ■ » 

Before proceeding further, equations 3-11, 3-12, and 3-13 will be 

written in finite difference form then rearranged in their final form as 

three linear, algebraic equations in the three unknowns o , a a , and v. Dur- 

r 0 

ing this process E and C are treated as variables which are functions of the 

location of c on the stress-strain curve. In the elastic regime, E and C 
eq 
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- i 

t 

are constants while in the plastic regime E is taken as the tangent 
modulus with C assuming the plastic wave speed based upon the current 
value of E. 

If we consider any two points 1 and 2 along a + c character- 
istic and apply equation 3-11 between them we get 

°rl * <* < v 2 • V ; 

• H [■ <D r * V + l> ” cv ] 1 + [- <°r ' °6 > + P ” CV ] 2 } 

where the right side of the equation has been averaged between points 
1 and 2 and the quantity c also denotes an average value between the 
same points . 

If we suppose now that the variables are known at point 1 while 
those at point 2 represent unknowns, we may transpose all variables with 
subscript 1 to the right side of the equation while retaining all sub- 
script 2 variables on the left side with the result that, 

. °r2 * pJV 2 • [ (o e • V r' + p >’ cv r/t r 2 ] 2 

. *L (a e -°T y ^r^ + pvcv ^r^J 1 + <, r l - pcv l 

If we clear parentheses on the left side of this equation, collect 

-terms in o , C. and v, and denote by a subscript 2 those variables which 
. r o 

t are to be evaluated at point 2, we obtain 
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f i + 1 

L r, *5 r„J 


A r 


a r2 


: 1 r 2 


- rpc 


+ PVC„ — 


Ar 


-1 


+ ^2^ £ 


' [<°e ' °r > + ‘ ,VC v ]j r'j + r 2 + °rl * p “ v l 


4-1 


Since the equation along a-c characteristic differs from (4-1) only 
in two algebraic signs, an identical development of equation 3-12 along 
a-c line leads to 

[ l ‘ r^+ rj a x2 + Tfh; a Q2 + [»' ‘ puc 2 v 2 

" * [ (o e ' CT r ) ' ( " ,cv ] 1 r 2 + °rl + ppv l 


4-2 


-In analogous fashion equation 3-13 along a dr = 0 line may be written 

AS 

it * '** (v i + v 2> i 

°r2 ' °rl " v L a 92 * V " ~ 2 J 


where E denotes the average value between two adjacent points. 

t 

By expanding the above equation as before, collecting terms in 

* 

a i o c , and v, and placing- subscript 1 and subscript 2 variables on 

r 6 

opposite sides of the equation we have, 


1 E At 

- T cr-„ + r— • v„ = a 


E At 


“x2 v 02 2r u ’2 rl v 02 2r u 1 


4-3 


Referring again to Figure (5), let (i,j) be a general point 

at which the variables o , a. , and v are to be solved. Then the ad- 

r 0 

jacent point along a + c characteristic at which the variables are 
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knowit-ts- -labeled (i-l-,-j)-i— Similarly, the adjacent point along a - c 

characteristic has indices (i+l,j-l), and finally along the dr = 0 line 

the adjacent point is denoted by (i,j-l). 

With this notation equations 4-1, 4-2, and 4-3 can be written 

"in their final form, with subscripts for a , a a , and v enclosed in 

t o 

parentheses. — 


[l + r & V - ~1 o (i ,j) a Q (i, 

L r i-l + r i J r r i-l + r i 9 


j) 


- + vc . 1 V (i, 

L 2 r t-1 + 


j) 




+ o r (i-l,j) - p 


+ c i,j 


v(i-l,j) (along + c) 

4-4 


[ Ar 1 . Ar 

’ J V 1 - 


Ar 

i + l i+l 


. r c i»3 + c i+l,j-l 

+ P L 2 


~ uc. 


Ar *1 

_ — J v(i,j) 

i i+1 J 


- [o e (i + l,3-D - o r (i+l,J-l) - pvc 1+1 ^ vci+lj-l)] - -f - r — 


+ a r (i+l,j-l) + p 


c i,j + c i+l,j-l 
v(i+l, j-1) 

(along 

4-5 
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1 E i,j +E i,j-1 

.(i.J) - At v(i,j) 


v ~Q 


4 r^ v 


E i,j + E i,j-1 

- V 1 ’*- 1 * - - - - - h ~ 7 


At v(i,j-l) (along dr = 0) 

4-6 


The foregoing relations constitute 3 equations in 3 unknowns 
and are of the form AX = B where [x] is the column matrix of unknowns, 

[A] is the coefficient matrix, and [B] is the column vector corresponding 
to the right side of the equations. The system as given by 4-4 through 
4-6 is valid for an interior point of the characteristic grid; however, 
along the i = 1 line a is given as the applied pressure vs. time. In 
this case only two equations along the — c and dr = 0 lines are used. 
Thus, o r is transposed to the right side of the equations as a known 
quantity and equations 4-5 and 4-6 are solved as 2 equations in the two 
unknowns and v. 

* 

Similarly, when i = NA at the outer edge of the disk, only the 

» 

equations along the + c and dr = 0 directions are used with = 0. 

Again we have 2 equations. (4-4 and 4-6) in 2 unknowns <j and v. 

0 

The equations used to compute the variation in shock amplitude as 
a function of radial distance r for the general case have already. been 
briefly discussed. The finite difference form of the equations follows. 

Referring to Figure 5, we may develop equation 3-11 between 
points (i,j) and (i-l,j) along the upstream edge of the shock in a man- 
ner analogous to that leading to 4-1 except that the values of c assume 
only the elastic value c q . 
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r. . 1 .. .. & r .. .. 

L *l-i~ r ~i -I ' JJ ‘ V l,J ' 

* , . . » 

‘ - P c o [i + v j— ^- ] v (i,j) 

+ p u c o v(i-l,j)] 


+ a r (i-l,j) - pc Q v(i-l,j) 


4-7 


Similarly, equation 3-25 may be written 


a r (i,j) + pc o v(i, j) = a r (i,j-l) + pc o v<i,j-l) 4-8 


Finally, equation 3-27 can be written 


o r (i,j) -|o e (i,j) =■ cr r (i,j-l) 


4-9 


Equations 4-7, 4-8, and 4-9 thus comprise the unloading shock equations. 


Computational Procedures 

An outline of the computational procedures will be given next 

• 

as an aid to understanding the program computer listings given in 
Appendix D. 

In this problem either a smoothly increasing and decreasing 
pressure pulse as a function of time is considered or a ramp input 
followed by a constant plateau. 

Referring again to Figure 5 and Appendix D, after the initial 
data is read into the computer the following computations are performed: 
(line numbers refer to the computer listing in Appendix D) 
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1) A differentiation subroutine is called to compute the 
value of E (tangent modulus) at selected points along 
the specimen stress-strain curve (line 0046). 

2) The variable plastic wave velocity C(I) and the constant 
elastic wave velocity CO are computed at points corre- 
sponding to the previously determined values of E (lines 
0048 and 0061). 

3) . Lines of constant R are determined using the input incre- 

ment £R (line 0064). 

4) The intercept time of ah elastic wave with each of the 
constant r lines is determined (line 0067). 

5) All variables along this line except C and E are set = 0. 
These are set equal to their elastic values (line 0069). 

6) From point i = 2, j * 1, the time of intercept of an elastic 

wave with the 1=1 line is computed (line 0092). This is 

% 

point 1=1, J = 2 . 

I 

7) The Aitken-Lagrangian interpolation routine (ATSM and ALI) 

' is called to interpolate the value of input pressure, a r> 

at this point (line 0094). 

8) Equations 4-5 and 4-6 are called with GAUSS to solve. for 

and v at this point (line 0104). 

9) The equivalent stress is computed and a check is made to 
see if yielding or unloading has begun (line 0106). 

10) New values of C and E are interpolated, for recomputing 

the average values specified in the equation coefficients 

(line 0115). Equations 4-5 and 4-6 are solved a second 

time for "improved" values of o and v. 

0 
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11 ) Using the final values of c > CT a and v, the elastic 

r 6 

strains are computed (line 0139). 

12 ) Beginning with line 0159 the iterative, calculations for 
the plastic strain increments are carried out prior to 
adding them to the elastic strains (line 0182). See 
Step 15 for a more thorough discussion. 

13 ) From point i = 1, j = 2 the intercept of a + c character- 
istic with the i = 2 line is computed . This point is the 
initial estimate of TA(2,2) (line 0198). Similarly from 
point i = 3 i j = 1, the intercept of a - c characteristic 
with the i = 2 line is computed. This is denoted by TANN 
(line 0202). In general TAPP find TANN will not coincide 
in the plastic regime in which case the earlier arrival is 
designated as the actual tine intercept for that particular 
i, j point. In this example i = 2, j = 2 . If TAPP > TANN 
an interpolation subroutine (TRANSL) is called in which the 
values of the variables at point i = 1, j = 2 are reduced 
linearly by the ratio (TAPP-TANN)/ [TAPP-TA( 2, 1) ]. In this 
way a new value of TAPP originating from a new point corre- 
sponding to the values of the interpolated variables more 
‘closely coincides with the value of TANN. If the stress 
field were entirely elastic in this region the foregoing 
procedure would be equivalent to translating the character- 
istic line segment defining the TAPP intercept parallel to 
Itself until TAPP coincided with TANN. Of course, if the 
field is completely elastic then TAPP and TANN automatically 


coincide 



If TANN > TAPP a similar procedure is carried out 

with the interpolation being performed now on the variables 

along the i = 3 line. Of course, in this particular case 

when j = 2 an interpolation cannot be carried out since it 

would involve variables on the j = 1 line and also a j = 0 

line which is non-existent. At all other points, however, 

for which j > 2 the interpolation is performed (line 0218). 

The process described in Step 13 is continued for each i 

line through i = NA. The procedure then returns to the 

i *= 1 line with j incremented by one and is repeated until 

all of the input values of pressure vs. time have been used 

or until the unloading shock occurs (line 0112).- 

Whenever the value of a at a given point becomes > SIGMAX, 

eq 

which is the largest plastic value of attained thus far 

along a particular i line, the material is continuing to 

yield. An initial estimate is made of the equivalent plastic 

» 

strain increment A£p by taking the difference between the 
current and preceding interpolated values of plastic strain 
after the latter has been increased by the elastic strain 
increments (line 0164). This value is then used in equations 
3-21 to 3-23 to compute estimates of the individual plastic 
strain increments (lines 0168 - 0170). These are used in 
equation 3-20 to compute an improved estimate of Atp. The 
iteration is repeated until two successive estimates of the 
strain increments differ by no more than TESTEP, an assigned 
test value.. • 

Several checks are made at each point after the final 
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17 ) 


computation of equivalent stress. If a > a , i.e., 

. eq o 

the material is in the' plastic range at a given value of 
i, and if a at a point (i, j) becomes less than the 
value at point (i, j-1) then unloading has occurred and 
an elastic unloading shock propagates radially outward 
from this point (lines 0108 - 0112 typically). Unloading 
is assumed to occur first along the 1=1 line 
[SIGEQ < SIGMAX(1)3» If it occurs along some other I line 
an error message is printed and the program terminated. 

If ISWITCH = 1 transfer to the unloading shock computations 
is bypassed . 

The time intercepts of this shock with each i line are com- 


puted (line 0353). 

18 ) At each intercept point a search is made for bounding values 

of time, TA (line 0354). - 

19 ) The problem variables are interpolated between these values 
of time and are stored as the downstream ' (earlier time) 
values along the shock. If the shock line does not cross 
any of the j -lines and therefore there are no bounding 
values, then the last previous values are stored along the 
•downstream edge of the shock (line 0379). 

20 ) Equations 4-7, 4-8, and 4-9 are called to compute the values 
of the variables upstream of the shock. The process is con- 
tinued to the i = NA-1 line (line 0405). 

21 ) At each i line a comparison is made between the values of 

O downstream and upstream of the shock. If the down- 
eq 

stream value is less than the upstream value, the shock 
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• I 

computations continue , otherwise the shock computations 
are terminated (line 0453) . 

22) At the i = NA line a r — 0. The jump in the value of v 

across the reflected shock is computed from equation 3“26 

(■line 0461) along with the accompanying values of <y and 

0 

the strains. The shock computations are then terminated. 

23) All previously computed variables for values of time 
greater than the shock intercept times are set equal to 
zero (line 0368). 

24) If the shock does not propagate to the NA line but is 
terminated as discussed in Step 21, the time intercepts of 
the extended shock line are computed (line 0473) . 

25) At each of these points a search is made for the bounding 
values of time and the problem variables are interpolated 
(line 0477) as in Steps 18 and 19 and are stored. All later 

. ^ i 

values previously comppted, if any, are set equal to 0. 

26) The calculations then return to the 1=1 line and the 
usual characteristic computations continue in the elastic 
mode until the values of input pressure are exhausted. The 
final stresses and strains thus computed are the residual 
Stresses and permanent set strains. 
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SECTION V. 

• ‘ » - 

EXPERIMENTAL EQUIPMENT 

As mentioned previously, a modest experimental program con- 
sisting of two methods for initiating a pressure pulse was conducted 
to verify the theoretically predicted solutions for the distribution 
of strain history in the specimens. 

A schematic diagram of the experimental fixture used in the 
hydraulic technique is shown in Figure 1 while Figure 6 presents a 
photograph of the fixture and instrumentation. The explosive method 
utilizes only the support plate and a tripod supported steel rod cover- 
ing the specimen central hole to retain and better distribute the ex- 
plosive gases. 

A pressure pulse is initiated on the inside surface of the 

% 

specimen by igniting a quantity of explosive with a fuse or by dropping 
a weight onto a plunger in the case of the hydraulic method. Since it 
has been stated that the hydraulic method produced poor results, the 
following equipment specifications and configurations are those used 
with the explosive method. 

Two 350 ohm strain gages are mounted on the surface of the disk 
—at a radial distance of 0.75 inches from the center. One gage is 
oriented with its sensitive axis along the radial direction while the 
second is mounted with its sensitive axis oriented perpendicular to the 
first in order to measure circumferential strain. The gages are connected 
to a bridge circuit and amplifier which in turn transfers the signals to 
a 2-channel oscilloscope where they are permanently recorded on 
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Polaroid film. 

The pressure pulse is sensed by a Piezo-electric pressure 
transducer and integral amplifier unit before proceeding to a single 


channel of a dual channel oscilloscope where it is also recorded on 

film. r . t : .... 


Listed below are the manufacturers, model or type numbers, 
and major performance parameters of the various instruments used with 
the explosives technique. 


Performance 

Item Manufacturer Model /Type Parameters 


Strain Micro-Measurements Type: 350 ohm 

Gages Company ED-DY-125AD-350 gage factor = 3.20 

Pressure Piezotronics Model: 109A pressure 0-80000 

irTransducer : . : ~ r : : psi; frequency 

range 0-500000 


hertz; rise time 
1 micro-sec. 


Bridge Leeds & Northrop * - — 

Company 

t 

Bridge Vishay Instrument --- 

Amplifier Company 

Oscilloscopes TEKTRONIX Model: 551 

dual beam 

. Plug-in Type: 53/54C 

Unit 


frequency range 
0-120,000 hertz 

frequency range 
0-120 ,000 hertz 


pre-amp rise time 
0.01 micro-sec. 



SECTION VI 


RESULTS 

* \ 

The hydraulic tests were the first ones conducted and re- 
sulted in typical pressure histories and strain histories as shown in 
Figures 7 and 8. As can be seen in Figure 8, the circumferential strain 
initially is compressive before becoming tensile which is contrary to 
expectations. Furthermore, the time at which the first positive pulse 
of the radial strain peaks is much too long relative to the peak of the 
pressure pulse. As a result there was some doubt as to whether the re- 
— corded strain signatures were due solely to the pressure pulse. It was 
anticipated that plate bending modes may also have been excited, es- 
pecially in view of the bending moment applied by the O-ring seals dur- 
ing clamping. 

As a result of these considerations it was decided to conduct 
the explosive test program to obtain more reliable data. The strain 
histories from the explosive tests presented in Figures 11 and 12 show 
much better behavior than those of the hydraulic method. Accordingly, 
the results obtained from the hydraulic method will not be discussed 
further while' those of the explosive program and two other test cases, 
as well as the multiple runs option on a prestressed specimen will be 
considered subsequently. 

Three test cases have been evaluated as a means of verifying 
the validity of the theoretical solution. In addition, a fourth case 
has been' run to illustrate the use of the multiple runs option in which 
a set of final permanent stress and strain distributions and equivalent 
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yield stresses for each I -line are used as initial values for a 
succeeding run. • ' 

The cases considered include: 

a) Steady state elastic solution with specimen 1 of 
Figure lb 

b) Steady state elastic-plastic solution with specimen 1 
of Figure lb 

c) Transient solution using the pressure pulse of Figures 
10 and 12 with specimen 2 of Figure lb 

d) Transient solution using a previously stressed specimen 

. (specimen 1 from Figure lb) with the pressure pulse 

. shown in Figure 25 . 

Each of these cases will be discussed individually. In each 
Instance, a positive stress denotes compression while a negative stress 
denotes tension. 

> 

Case a) The steady elastic solution has been obtained by applying a 
pressure input consisting of a ramp to 5000 psi in 8 p,sec followed by 
a constant plateau of 5000 psi to 225 ysec and program termination. 

A fictitious yield stress = 20000 psi has been used to insure an elastic 
stress field. The final stresses at four selected radial locations have 
been obtained from Figures 13, 14, 15, and 16 and are tabulated below 
In addition to the values calculated from the standard thick walled 
cylinder equations from reference 12.^ It is noteworthy that during 
the first 25 to 35 (ji seconds, depending upon the particular I-line, 

^The thick walled cylinder equations referenced in cases (a) 
and (b) assume a state of plane strain rather than plane stress. These 
equations nevertheless serve as a useful approximate check. 
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these curves behave in the same manner as those given in reference 8. 
Thereafter the disturbing effects of the reflected waves from the free 
boundary are evident in the form of oscillations about the estimated, 
best-fit steady state solution. These appear to continue indefinitely 
in the absence of damping effects. 



1*1 

1 = 3 

1 = 7 

I = 19 


R = 0.75 

R = 1.00 

R = 1.50 

R = 3.00 

Computer solution a r : 

5000 

2800* 

1200* 

200* 

a e : 

-5000* 

-2800* 

-1400* 

-500* 

Thick Walled Cylinder a v '■ 

5000 

2750 

1140 

178 

• c e : 

-5260 

-3020 

-1425 

-464 


^Estimated best fit steady state solution. 

The agreement between the computer solutions and the known 

solutions is quite good, on the order of 5 to 10 percent. 

% 

Case b) The steady state elastic-plastic solution has been obtained in a 
manner identical to that of the previous case except that the actual 
yield stress of 2500 psi has been used. 

In this case the circumferential stress does not attain as 
* large a negative (tension) value as in the elastic case. Furthermore, 
because of the slower plastic wave speed the effects of the reflected 
vaves do not arrive at the 1=1 line until a later time. 

The final stresses at the same radial locations as case (a) 
have been obtained from Figures 17, 18, 19, and 20. They are tabulated 
below in addition to the values calculated from the elastic-linear 
. strain hardening, thick walled cylinder equations given in reference 11 
using an equivalent yield stress of 4000 psi. Again agreement between 
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the computer solutions and the closed form solutions is reasonably good 
with an error on the order of 5 to 10 percent except for the circum- 
ferential stress at the 1=3 line. Also note that it is somewhat 

difficult to accurately estimate the best-fit steady state response in 

2 

the oscillating region of the curves. 



I = 1 

1 = 3 

1 = 7 

I = 19 


R = 0.75 

R = 1.00 

R = 1.50 

R - 3.00 

Computer solution a : 

5000 

3600* 

1800* 

300* 

a 0 1 

-500* 

-800* 

-2000* 

-800* 

Thick Walled Cylinder a : 

5000 

3740 

2060 . 

280 

; 

-415 

-515 

-1960 

-800- 


♦Estimated best fit steady state solution. 

• . 

Case c) Figures 21 and 22 present the stress-time computer solutions 
along the 1=1 and 1 = 24 lines using the typical pressure pulse of 
Figures 10 and 12 while Figures 23 .and 24 present the strain -time 
solutions for the same data at the same radial locations .' The I = 24 

f 

line (R = 0.75 in.) corresponds to the location of the strain gages on 
the experimental specimen. These provide the only means of correlating 
the experimental and theoretical results for this test data. Figure 11 
presents typical tracings of the experimental strains -while Figure 12 
shows a plot of the same data. The plotted data has also been super- 
imposed on the computed solutions of Figure 24 for ease of comparison. 

The information required to convert the experimental pressure 
and strain tracings to actual pressure and strain values is given below. 
The calibration factors were taken from the appropriate instrument 

2 

See Footnote (1) page 
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manufacturer’s instruction manual. 

Oscilloscope Settings During Explosive Tests, 



Vertical Sensitivity 

Sweep Rate 

- 

(mv/ cm) 

(psec/cm) 

Pressure transducer 

200 

2 

Radial stress gage 

50 

2 

Circumferential strain gage 

50 

2 


The slope of the pressure input/voltage output calibration 
curve of the pressure transducer is 18.3 psi/mv. 

A known test signal from the strain gage bridge instrument was 
used in conjunction with the manufacturer's calibration factor to con- 
vert the strain gage data to actual strain values. This conversion 
factor is 1.95 p in/in per mv. 

Comparison of the experimental data in Figure 24 with the com- 

% 

puted solutions indicates that the magnitude of the first positive peak 
of the measured radial strain is 56 p in/in while the computed value is 
42 p in/in. The corresponding values for the second positive peak are 
51 p in/in versus 15 p in/in respectively. It can be seen that the 
times at which the measured and theoretical peaks occur agree exactly. 
Similarly, the first negative peak of the experimental results agrees 
precisely with the computed value, whereas the second negative experi- 
mental peak is -55 p in/in compared to the smaller theoretical result 
of -15 p in/in and occurs approximately 1 psec later. 

The behavior of the circumferential- strain is predicted with 
reasonable accuracy; the major difference being in the magnitude. 
Experimentally, -28 p in/in was obtained versus -18 p in/in for the 
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predicted value. The. times at which these peaks occur are again seen 
to be identical. * 

Although there are some discrepancies between the experimental 
and theoretical results as discussed above, the shape and phasing of 
the curves is excellent. In view of the somewhat unsophisticated ex- 
perimental technique used it is also considered that the magnitude of 
the curves are in reasonable agreement. In particular, the discharge 
of a quantity of powder in the central specimen hole creates a pressure 
pulse that is not entirely uniform over the depth of the hole. As a 
result, a purely radial pressure pulse is not realized, but instead a 
complicated system of oblique waves resulting in reflections from the 
—specimen free surfaces, the generation of shear waves, and complicated 
wave interaction effects. 

In view of these considerations it is concluded that the ex- 
perimental data presented is in agreement with the theoretical calcu- 
lations to within the expected accuracy of the experimental technique. 

Case d) This run has been conducted to illustrate the analysis of a 
previously stressed specimen wherein each I -line begins with permanent 
set values of the stresses, strains, and yield stress obtained from a 
previous run. These are read in either as initial values along the 
J “ 1 line or merely as values associated with each particular I-line. 

To illustrate, the following arbitrary but typical values of the 
variables were input as initial data for specimen 1 which was subdivided 
into 32 I-lines. 
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Variable 

1-line Range 

Value 

EPSEQl(I) 

1-7 

.0021265 


8-32 

‘ .0002500 

SIGMAO(I) 

1- 7 

4823.8 

■ ^ 

8-32 

2500.0 

STNRA(I.l) 

All 

.0001 

STNCA(I.l) 

All 

- .0001 

STNZA (1,1) 

All 

- .0001 

SRADA (1,1) 

All 

500 - 

SCIRCA(I,1) 

All 

” - 500 


Typical results from this run using a 7000 psi input pulse 
are shown in Figures 25 and 26 for the stress and strain distributions 
along the 1=1 line. , 

As can be seen from the plots, the initial values of the 
variables are those input as the initial ‘conditions . Thereafter the 
stress and strain time histories proceed in a rather normal manner but 
are displaced positively or negatively in accordance with the initial 
permanent set values. • 
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SECTION VII 

CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER WORK 

Section VI has presented the results of the computer studies 
which have been run to verify the validity of the theoretical analyses 
in comparison with known solutions and experimental test results. 

The elastic steady state and elastic-plastic steady state solu- 
tions agree quite well with known closed-form solutions. In addition, the 
test results agree with the theoretically predicted solutions to within 
the expected accuracy of the experimental techniques. The case showing 
the multiple runs option with a prestressed specimen yields results that 
would be expected. 

In view of the above, the following recommendations for addi- 
tional work are suggested: 

1) Extend the program to include other similar geometries such 
as spheres and cylinders with various end conditions . 

2) Incorporate the capability to evaluate strain rate effects. 

- 3) Add a subprogram to predict reverse yielding during unloading 
and to compute the associated stress and strain distributions. 

In final summary, it is concluded that the method of analysis 
presented in this work and the accompanying computer programs provide a 
means of calculating the transient stress wave and unloading shock wave 
behavior in a disk or open-ended cylinder in which the plane stress 
assumption is not greatly in error. Of even more importance, it lends 
some validity to the assumption that this general approach could be 
applied to any problem involving propagation of one-dimensional stress 
waves under states of combined stress in the elastic-plastic range. 
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SELECTED NOMENCLATURE AND COMPUTER NAMES 

Symbol Computer Name Description 

A, B. Points on r-t shock diagram (Figure 4). 

Coefficient matrices. 


(bar) 

c CO, CA 

D - 

6 < ) 


Ac p « de p 





1 

J 


DEPP 

EO, EA 

STNRA 

STNCA 

STNZA 

I 

J 


V 

P 

r 


P 




RA 

RHO 


Denotes average of two adjacent values. 

Elastic or plastic wave speed. 

Point on r-t shock diagram (Figure 4). 

Denotes jump or discontinuity in quantity 
in parentheses. 

Equivalent plastic strain increment. 

Plastic strain increment. 

Modulus of elasticity or plastic tangent 
modulus . 

Radial strain. 

f 

# 4 

Circumferential strain. 

» 

Axial or longitudinal strain. 

Radial coordinate line number. 

+ c Characteristic number. 

Constant relating plastic strain incre- 
ment and deviatoric stress tensor. 

Poisson's ratio. 

Pressure. 

Radial coordinate distance. 

Mass density. • 

Deviatoric stress tensor. 

Radial stress. 


SRADA 



Symbol 



t 

u 

v 



Computer Name Description 

SCI RCA Circumferential stress. . 

— Axial or longitudinal stress. 


SIGEQ 


S1GMA0 

XA 

VA 


Equivalent stress. 

Mean stress. 

Yield stress. 

Time. 

Radial displacement. 

Particle velocity. 

Initial elemental volume. 

Final elemental volume. 

Axial or longitudinal displacement. 

Subscripts denoting adjacent character- 
istic points or variables on either side 
of an unloading shock. 
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ADDITIONAL SELECTED COMPUTER NAMES 


Computer Name 
All 
ATSM 

CINTRP 
DE1EP 
DE1RA 
DEPRP 
DEPCP 
DEPZP 
DERIVL 

EINTRP 
EPSEQ 
EQSEQI 
EPSIIN 

GAUSS 

* ISTART 
ISWST 

ITERM 

N 

NA 

NE 

NT 


Description 

Aitken-Lagrangian Interpolation Routine. 

A subroutine which orders tabular data suitable 
for use in ALI . 

Interpolated wave speed c. 

Increment between SIGMAE and EPSILN entries. 

Increment in radial coordinate r. . 

Radial plastic strain increment. 

Circumferential plastic strain increment. 

Axial or longitudinal plastic strain increment. 

A subroutine which computes the derivatives of a 
tabulated function. ' 

Interpolated elastic or plastic modulus. 

> 

Equivalent plastic strain. 

I 

Initial estimate of equivalent plastic strain. 

A tabular entry of strain from the stress-strain 
curve . 

A subroutine which solves systems of linear equations 
by Gaussian Elimination. " 

I-line at which unloading shock starts. 

Internal indicator controlling when strains during 
unloading are to be computed. 

I-l,ine at which unloading shock terminates. 

Integer index. 

Number of radial coordinate lines in characteristic 
-grid . 

. Number of entries of tabulated stress-strain curve. 
Number of entries of tabulated stress-time input. 



Computer Name 

RAl 

SIGE 

SIGMAE 

SIGMAT 

SIGMAX(I) 

SIGN 

SIGT 

TA 

TANN . 

TAPP 

TESTEP 

TIM 

TRANSL 

XEPRP 

XEPCP 
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Description 

Inside radius of disk. 

Stress in called subroutine variable list. 

Tabulated stress values from stress-strain curve. 

Tabulated stress values from stress-time input. 

Maximum equivalent stress attained along the ith r 
line . 

A subroutine which assigns the algebraic sign of one 
variable to another variable. 

Stress in called subroutine variable list. 

Time. 

.The trial time intercept of a - c characteristic 
'segment on an r line from the r + 1 line. 

• 

The trial time intercept point of a + c character- 
istic line segment on an r line from an r - 1 line. 

A small test number used to determine the convergence 
of an iterated variable. 

Time in called subroutine variable list. 

A linear interpol’ation routine used to effectively 
translate the variables associated with the larger 
'of the. TAPP or TANN intercepts to coincide with 
those of the smaller. 

Intermediate value of radial plastic strain incre- 
ment during iteration. 

Intermediate value of circumferential plastic strain 
increment during iteration. 

Intermediate value of axial or longitudinal plastic 
'strain increment during iteration. 


XEPZP 








Forces 
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Fig. 5 . --Characteristic Network with 
Indexing Notations and Applied Pressure Schematic. 
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I 



Fig. 6 . --Experimental Fixture and Instrumentation 

(Photograph) 
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i-g . 7. --Pressure Pulse Experimental 
l Specimen 1 of Figure lb (Hydraulic Method). 
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APPENDIX B 


CHARACTERISTIC DERIVATIONS 


From 

Section 

III equations 3-1, 

3-5, 

and 3-6 are 

- 


ha r . °r • "e 
dr + r " 

dv 

P 

- 



•S 1 

&v » i r i _ ■ 

Ba e'i 




dr E Ldt v 

dt J 




v _ i r^fa 




- 

r E. Ldt " v 

5Tj 

• 

Writing this 

system in determinant form along with the 

differential changes 

in a , a fi , and v gives 
r o . 


* a r 



- 

dv • _ c 0 " CT r 

Sr 



• 

p dt r 


i K 


E A 


E dt 

" E dt ’ dr 

“ U 

• 

v ScT r 

+ iiz§. 


V 


“ E dt 

E dt 


r . 

♦ 





3a * da 

r , r 

, _rr-_ dr + -r— 
dr dt 

dt 


“ da^ 



do 3c fi 

-if dr + ^f dt 


do e 


“• dr + ^7 dt = dv 
dr dt 
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Sa r 

Solving for — by Cramer's Rule yields 
of 


vhere 


B 


d<T r 

"dr 


A 

B 


1 

0 

0 

dr 


- 1 - 

E 

v 

E 

dt 


0 

0 

0 

dr 


v 

E 

l 

E 


dt 


0 ‘ p 

- 1 •: o 

o o 

0 V5 o 

o o 


dr 


dt 


and A is given on page . Expanding B on Column 3 gives 


B s dr 


1 

0 

0 

dr 

0 


1 

E 

v 

E 

dt 


JL‘ 

E 

1 

E 


P 


1 

0 

0 

dr 


P 

0 

0 

0 

dt 


Continued expansion on Column 1 provides the following re- 


duction 
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By expanding the first on row 3 and the second on row 1 we arrive at 


B - (dr)(l)(dt) 


1 10 

- 

i - 2. . ! 

E E 

i 0 0 

+ (dr) 2 (- p)‘ 

- H I o 

E E ■ 

0 dr ‘ dt 

- 

0 0 dr 


or finally, 

B - (dr)(dt).[+ |dt] 


x-) 2 Kt) - - 5 -] 


Clearing of brackets and setting the result to zero gives 

1 '2 3 /IN 2 v 2 3 

B - + ~ d r(dt) Z - p(dr)- 3 (jj + p — (dr) J = 0 . 

E 

* If we factor the common term dr we obtain 


dr[+ | dt J + dr 2 (- rj + 0 rj)] - 


L . E \ E 2 K E 

This equation is satisfied only when 


dr = 0 





The second relationship can be written 





therefore 



dr 

dt 


±r ^ 

p(l - v 2 ) 


•] 


% 


and dr = 0 . 


A «= 


Similarly, determinant A is 



Expanding on Column 3 
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Expanding the first on Column 4 and the second on row 4 yields 
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Clearing parentheses and setting A - 0 gives 


2 Cq - CT r 


*= dr dt — + p dr dv dt - - p dr da —r 

• r E Y E r p 2 


' 2 

. j 2 'v •*. v ,2 v , 
-k.p : dr - dt - + p dr — da r 


Fdt 2 — — + p dv dt — p dr da “ + p dr y dt — + p dr da ”] - 

EL- r r E it ' u r j 


either dr = 0 or the bracket =0. In the latter case, rewriting 


o fi - a _ da 2 

— dt Z + p dv dt = p dr — — - p dr v — dt - p dr da 

r K H E K r K E r 


Multiplying by dr and divide by dt 


<°e ' V dr . . dr /dr \ 2 r 

r ■ + p dv "dt ~ VdEV L 


2 r da 2 

[p -f - P ^ - p J- daj 


since ™ = ± c _ or 


* 2 


09 ■- 


2 . E _ 2 f . 2 V 

c and — = c (1 - v ) 


<°<S - °r> , , , 2 r da r v „ ” 2 do r *1 

— dr ± pc dv - c l~ 5 - -.P » 7 ^ 57 J 

c (1 - V ) c (1 - V ) 


[ 1 — (l - p2 ) - p » 7 dt ] 

L c (1 - v ) r J 
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< a e ~ a r > r 

V 1 — A*- j. A.. = I - 


X MM /}«* CS I *» •< — /» ^ J **”] 

vi — y w» v L. W ^r p 1/ ^ V UW j 


Replacing dt by ± — gives 
« c 


<o-q “ a r ) “ ± pc dv = [da r + puc — drj 


or finally, rearranging 


da r + pc dv = <Cg 


- a 


± pucv) ~ 


3-11 and 
3-12 


The top sign refers to propagation along the + c div action and the 
bottom sign to propagation along the - c direction. 

Along dr - 0, there is no wave propagation therefore we can ex- 
press a static relation that holds between the stresses and velocity 
from equation 3-6. 


v 

r 


E “ dt 
r 


v v dt 
E r d0 


i do e 


da 

v r 


E dt E dt 

» 

da. - v da 

v7 X 


= 1 - u 


da^ 

da- 


da. 

da 


£-[l - E “ “ 1 — 
£> L r de J v 


or 


da, 


r c v [ da Q - E ? dt ] 


3-13 
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APPENDIX C 


COMPUTER PROGRAM INPUT 


Following is a list of the program input, the card format, 
and a definition of each item with clarifying discussion as required. 


• Format (8F10.0) 

DELRA - radial increment of disk (in.). This quantity controls the 
characteristic mesh spacing and also strongly controls how 
...... rapidly the solution progresses. Too large a value could re- 
sult in unacceptable accuracy of the finite difference equations 
; . while too small a value will cause excessively small increments 
: - of the input pressure to be interpolated thereby increasing 

machine running time. In the present studies values from 1/8 
to 1/64 in. have been used. 

RAI - inside radius of disk (in.). 

7 ribm sec 3 ! 

RHO - specimen density — 7 - j J 

E0 - modulus of elasticity (lbf/in 3 ) 

NU - Poisson’s Ratio 

DEI£P - the increment between tabulated values of EPSILN(I). In the 

differentiation subroutine used in this work (DERIVL) this 
increment must be a constant. 

TESTEP- a test value used to test convergence of the plastic strain 
increments during iteration. In the present work .00002 was 
used . 


ei 


Format (615) 

NA - - 1 + number of radial increments into which the disk has been 
divided. It also equals the number of I lines. 

NE - number of tabulated entries in stress-strain arrays 
CSIGMAE(I) - EPSILN(I)]. 

NT - number of tabulated entries in stress-time input 
[SIGMAT(I) - TIME (I)]. 

NRUNS -^an indicator- showing whether the current run is the first with 

unstressed material or a succeeding run with prestressed material 

from a previous run. For the former case set NRUNS = 1; for the 

latter NRUNS = 2. See also discussion of SIGMA.O(I). 

JUM - an indicator which terminates the program when J £ JLIM. Used 

either as a normal termination condition or to prevent "runaway". 

KK - an indicator which controls whether or not plastic wave speeds 

will be computed. If KK = 1 only an elastic analysis is con- 

% 

templated and plastic moduli and wave speeds are not computed. 

» 

If KK 1 plastic moduli and both elastic and plastic wave speeds 
are calculated. -• 

INPUT - an indicator used to bypass the unloading shock calculations, for 
example, when the input consists of a "ramp" followed by a con- 
stant. plateau and one wishes to ensure that inadvertent round-off 
or interpolating errors do not cause SIGEQ to be less than 
SIGMAX(I) and thereby cause unloading. If INPUT = 1 the shock 
calculations are bypassed. If INPUT ^ 1 the shock calculations 
will be allowed to occur. 
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SIGMAO(I) - the stress at the proportional limit where the stress- 

I - ; 

strain curve ceases to be linear. This value is read in as an 
'array; one for each I line. When an analysis of a previously 
- unstressed disk is to be done each SIGMA.O(I) is assigned the 

value of the proportional limit. However, when the disk has 
' been previously stressed into the plastic range, some or all 
of- the I-lines will have experienced stresses above the pro- 
- r portional limit. The maximum stresses for each I-line are 

printed out as SIGMAX(I). Those greater than the original pro- 
portional limit should be read in as the values of SIGMA.0(I) 
prior to a succeeding run. They become, in effect, new pro- 
portional limits for their respective I-lines. When SIGMAX(I) 
is less than the original proportional limit, then the SIGMAO(I) 
for that line should be the original proportional limit. 

STNRA(I), STNCA(I), STNZA(I) - these are the residual permanent set 

* 

strains and are read in as 0 prior to an initial problem run. 

t 

For succeeding runs the final values from the previous run 
which occurred after the disk was completely unloaded are to be 
used. 

EPSEQI(I) - equivalent plastic strain initially computed by the program 
as the strain corresponding to the proportional limit. The 
final values corresponding to SIGMAX(I) are printed as output 
-and are to-be read in as initial values for succeeding runs. 

Format (8F10.0) 

TIME (I), SIGMAT(I) - tabular values of input time and corresponding 
pressure stress on the inner edge of the disk. 
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EPSILN(l)> SIGMAE(I) - tabular values of strain (EPSI1M) and corre- 
sponding stress from the stress-strain curve of the material. 
Only values for the plastic range are to be read, i.e., above 
the proportional limit. 

• . Format (815) 

IBRANCH - an indicator controlling the input of initial value data. 

If set — 1, initial values of the yield point, 3 strains, the 
Strain corresponding to the yield point, and 2 stress components 
are read in for each I-line. If < 0 the yield point is set to 
2500 psi with all other values =0. 

The following input data is listed at the end of the program 
and is used when output plots of the stresses and strains as a function 
of time for selected I-lines are desired. 

Format (615) 

NPL0T1 - the number of I-lines for which plots are desired. 

t 

IPLOT(Il) - the specific I-line numbers for which plots are desired. 

Format (8 f10 .0) 

XMIN - the minimum value of time utilized in the plot. This would 
usually be 0 . 

BELT - the increment in time per inch of plot space available. 

YMINl - the minimum values to be plotted of SRADA, SCIRCA, and VA on 
the same graph. If positive and negative values occur then 
YMINl should be chosen to at least encompass the largest nega- 
• 7 - tive number. 

DELYl - the increment in YMINl per inch of plot space available. 
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YMIN2- - the minimum values to be plotted of STNRA, STNCA, and STNZA 

on the same graph. If positive and negative values occur then 
YMIN2 should be chosen to at least encompass the largest nega- 
tive number. 

DELY2 - the increment in YMIN2 per inch of plot space available. 



APPENDIX D 


COMPUTER PROGRAM LISTING 


The computer program listing and sample output pages are not 
included in this report. 
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The natural frequencies and corresponding mode shapes for composite, 
axisymmetric bodies of varying area undergoing longitudinal vibrations 
are determined both experimentally and theoretically. Composite rods, 
cones and conical shells of Lucite, Nylon, Polycarbonate and Polyethylene 
are investigated. Strain measurements are taken from strain gages 
mounted on the surface of the models which are excited at varying fre- 
quencies by a shaker table. Lumped parameter approximations of the 
continuous bodies are solved numerically on the digital computer for both 
fixed-free and free-free boundaries in order to determine the natural 
frequencies, mode shapes and relative strains. The method of character- 
istics is utilized in a wave propagation approach to obtain the response 
of a point in the model to a random input. The random input selected 
is white noise. The natural frequencies are obtained by performing a 
Fast Fourier Transform Analysis on the response. 


LIST OF SYMBOLS 


INTRODUCTION 


A cross-sectional area 

a Lagrangian coordinate 

a.. influence coefficients 

»J 

c shift rate 

E modulus of elasticity 

fj fundamental frequency 

9 gravitational constant 

K spring constant 

L length 

t time 

u particle velocity 

x distance from base 

e engineering strain 

p density 

o engineering stress 

tp impact function 

* Graduate Assistant 
t Associate Professor 
£ Associate Professor 


In recent years engineers have been re- 
quired to design ballistic range models which 
are capable of being launched at hypersonic 
velocities for aerodynamic tests of ballistic 
missiles and space craft. At such velocities 
the designer is confronted with the dual prob- 
lem of using materials which are structurally 
strong enough to withstand the severe loading 
conditions and materials which minimize the 
ablation of the missile nose due to aero- 
dynamic heating. These models are frequently 
used to study boundary layer transition which 
could be influenced by modal vibration. The 
purpose of this paper is to study the vi- 
bration charac teri s t ics of the model as it is 
being tested in the range, and to develop an 
analytical means of predicting the transient 
response of the model under dynamic loadings. 

The problem of vibration in continuous 
systems such as cylindrical rods has been in- 
vestigated extensively in various research 
papers and textbooks £l-4]. However, vehicles 
are generally made up of conical shells, sec- 
tions with abrupt changes in area, and lam- 
inated sections of different materials. There 
is a definite lack of experimental and 
analytical developments concerning vibration 
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and wave propagation in structures of varying 
area and/or multi-layered with changing im- 
pedances. 

Experimental tests were performed using a 
shaker table to excite the models. Models 
were made from Plexi-glass (Lucite), Nylon, 
Polycarbonate (Lexan) , Polyethylene, and the 
combination thereof. Strain gages and accel- 
erometers were mounted on various positions on 
the model to determine the natural frequency, 
mode shape, and transient responses in the 
model. Two theoretical approaches, classical 
and wave propagation, were utilized. Natural 
modes of vibration and mode shapes were deter- 
mined by first the classical approach, which 
constitutes the solution to the boundary value 
problem and in turn the eigenvalue problem. 
Numerical solutions were obtained by using a 
system of spring-masses. In the wave propa- 
gation approach the method of characteristics 
was utilized. A numerical integration of the 
characteristic equations for the range model, in 
conjunction with appropriate procedures to en- 
"sure that the conditions on the external boun- 
daries as well as the continuity conditions at 
the sections with the step changes in area and 
varying impedances are satisfied, will yield 
-the stress, displacement, and the particle 
veloci ty anywhere in the model Q>J. Since this 
method is applicable for any type of loading, 
the natural frequencies can be obtained by 
applying a random input to the model and then 
performing a Fourier analysis on the response 
of some point in the model. 

EXPERIMENTAL INVESTIGATIONS 

A variety of models including rods, com- 
posite rods, composite cones and composite con- 
ical shells were fabricated and tested. 

Sketches of these models are shown in Fig. 1. 
Plastic materials were utilized rather than 
metals in order to maintain lower natural fre- 
quencies. An upper limit on the frequency of 
the forcing function was imposed by the shaker 
which was used - MB Electronics Model PM-50. 

All models were constructed from solid 
plastic rods which were turned on a lathe to 
their final dimensions. In the case of a com- 
posite model, the sections were glued together 
using an epoxy adhesive. Strain gages were 
then mounted on the surface and the model was 
glued to a base plate which was bolted to the 
shaker table. As the input frequency to the 
shaker was varied, the output from the strain 
gages was noted. In this manner the natural 
frequency of each model was determined. 

Since the elastic modulus of many plastic 
materials varies due to manufacturing incon- 
sistencies, age and environment in addition to 
the effects of turning on the later, exper- 
-iments were performed in order to determine 
the proper value of the modulus for the materi- 
als being used. Solid rods were used as the 
test specimen since the fundamental frequency 


equation for longitudinal vibration of free- 
fixed rods is well established. 



The fundamental natural frequency of each rod 
was determined experimentally as was the density 
Using this information and Eq. (1) the value of 
the modulus for each material was calculated. 
These experimentally determined values are 
given in Table 1, and are used in all subse- 
quent calculations. 

TABLE 1 


Experimentally Determined Values of the 
Modulus of Elasticity and Density 


Material 

E (lb/in^) 

p( lb/i n S ) 

Polyethylene 

2.8x]0 5 

0.0348 

Lucite 

3.3xlO S 

0.0420 

Nylon 

6.0xl0 S 

0.0430 

Lexan 

3.9xlO S 

0.0447 


A series of tests involving composite rods 
made up of two rods glued end to end were per- 
formed. These experimental results, given in 
Table 2 and Fig. 2, are restricted to the 
first two natural frequencies due to the upper 
frequency limit of the shaker. The results of 
the lumped parameter computer analysis which is 
explained in the next section is also given 
there. 

The next series of tests involved com- 
posite cones made up of a base and a tip glued 
together. Experiments were performed with the 
base and tip of the same materials as well as 
different materials. Table 3 and Fig. 3 
illustrate the results of these experiments and 
also the lumped parameter computer results. 

This also gave a check on the influence of the 
epoxy bond between the sections. Based on the 
results, it was concluded that the influence 
of the epoxy was negligible. 

The final experiments were performed using 
a concial shell as the base and a conical solid 
of a different material as the tip. These re- 
sults are shown in Table It and Fig. 4. 

LUMPED PARAMETER COMPUTER ANALYSIS 

The physical models are approximated by a 
lumped parameter system consisting of con- 
centrated mass elements and concentrated spring 
elements. Each model is divided into equal 
length segments with the mass treated as con- 
centrated at the mid-point of the segment. The 
rods are divided into 16 segments and the cones 
into 24 segments. The mass of each segment 
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varies as its volume in the case of the cones. 
The mass elements are connected by spring 
elements with spring constants determined from 
the properties of the material according to the 
equation: 


1 



* *(i) 

x(i-l) A(x 


( 2 ) 


where x is measured from the base and: 


K. j = spring constant between i th and 
(i-1) th positions 

x(i) = distance at i th position 
x(i-l) = distance at (i-1) th position 


The computer program utilizes the influence co- 
efficients, a . j , which are determined by 


’ij 


•j,-3 

J L = 1 


1 

K(L) 


(3) 


The computer program which is used to 
analyze the lumped parameter models first deter 
mines the eigenvalues by (1) pre-mul t ipl y ing 
the mass matrix by the influence coefficients 
to obtain the dynamic matrix, (2) transforming 
the dynamic matrix to a Hessenberg matrix. The 
natural frequencies are then determined from 
the eigenvalues. The eigenvectors and mode 
shapes are then calculated 

Natural frequencies obtained by this 
method are compared to the experimental results 
in Tables 2,3, and 4 and in Figures 2,3, and 4. 
The first and second mode shapes and corres- 
ponding material frequencies, for additional 
Composite cones, both solids and shells, are 
computed by this lumped parameter appraoch. 
These results are given in Table 5 and in 
Figures 5 and 6. The relative values of strain 
for the first and second modes are also com- 
puted and are shown in Figures 7 and 8. 


Table 2 


Composite Rods, Experimental vs. Lumped Parameter Frequencies 
(See Fig. 2 for dimensions) 


Material 

Frequency (KHz) 

Base 

Tip 

Experimental 

Lumped Parameter 



1st 

2nd 

1st 

2nd 

Polyethylene 

Luc i te 

1.6 

5.35 

1.64 

5.27 

Nylon 

Lexan 

2.0 

6.0 

2.16 

5.98 

Nylon 

Polyethylene 

2.25 

5.6 

2.28 

5.68 


Table 3 

Composite Cones, Experimental vs. Lumped Parameter Frequencies 
(See Fig. 3 for dimensions) 


Material 

Frequency (KHz) 

Base 

Tip 

Experimental 

Lumped Parameter 



1st 

2nd 

1st 

2nd 

Polyethylene 

Nylon 

4.45 

.io;2 

4.52 

10.17 

Polyethylene 

Luci te 

4.5 

— 

4.53 

9.21 

Polyethylene 

Polyethylene 

4.5 

9.2 

4.64 

9.26 
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5 ASSESSES 


Table 4 


Composite Cone, Shell Base with Solid Tip, Experimental vs Lumped Parameter 
Frequencies. 

(See Fig.' 4 for dimensions) 


Material 

F requency . (KHz) 

Base Tip 

Experi 

1st 

menta 1 
2nd 

Lumped Pa 
1st 

rameter 

2nd 

Polyethylene Lucite 

3.15 

8. 1 

00 

f*N 

J- 

oo 


Table 5 

Composite Cones, Lumped Parameter Frequencies 
(Base 1 in. dia. x 3 in. long. Tip 3 in. long) 


Materia 1 

Freauencv (KHz! 

Graph Symbol 

Base 

Tip 

1st 

2nd 


Poly, solid 

Poly, solid 

4.64 

9-25 

& 

Poly, solid 

Lucite solid 

4.45 

9-21 

■J- 

0.1 in. Poly 
shel 1 

Luci te sol id 

3.60 

10.06 

• 

0. I in. Poly 
shel 1 

Poly, solid 

3.80 

10.22 

o 

0.2 in. Poly 
shel 1 



Poly. soli<J 

4.38 

11.01 

□ 



Fig. 1 - Experimental models 


Fig. 2 - Composite rods, experimental vs. 
lumped parameter frequencies 
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lugnitud*) STRAW {r*l» siAgnltud*) 







STRAIN/MAX. STRAIN 



o .5 1.0 


DISTANCE 

Fig. 7 - Composite cones - first mode strain 
(see Table 5 for symbols) 


WAVE PROPAGATION ANALYSIS 

The stress wave propagation theory used in 
this analysis was developed independently by T. 
Von Karman CQ and G. I. Taylor [7] in 1942 
and used by many investigators in various forms 
since that time. It is a one dimensional 
elastic-plastic theory used to determine the 
transient properties of solid materials under 
dynamic loading. 

The response of materials is governed by 
the equations of continuity and momentum which 
are given in the Lagrangian coordinate system 
as follows: 


** 

Is 


W 

&a 

at 


In. c 

1 

3(<tA) 

(5) 

at 

PA 

aa 


From the above equations the characteristic 
equations may be obtained: 



where 

2 I dor 

C ' f JE (7) 

and 

dip = cde (8) 
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STRAIH/HUCIMUN STRAIN 



DISTANCE 

Fig. 8 - Composite cones - second mode 
strain 

(see Table 5 for symbols) 


A numerical integration of the characteristic 
equations in the range model, in conjunction 
with appropriate procedures to ensure that the 
conditions on the externa) boundaries as well 
as the continuity conditions at the sections 
with the step changes in area and varying im- 
pedances are satisfied, will yield the stress, 
displacement, and the particle velocity any- 
where in the model £5] • Since this method is 
applicable for any type of loading a random in- 
put was applied to the model. The random in- 
put signal selected is white noise which has 
a delta function as its autocorrelation 
function. The transient response at a point 
on the model was obtained numerically on the 
IBM 7094 computer. A Fast Fourier Transform 
method was performed on the response to obtain 
the natural frequencies of the structure. The 
advantages of. the Fast Fourier Transform are 
two-fold.- The number of actual arithmetic 
operations is reduced drastically, causing 
increases in speed of several orders of 
magnitude for reasonable record lengths. Also, 
because of the fewer operations performed, 
truncation and roundoff errors are reduced, 
producing a more accurate result. The results 
of the wave propagation analysis are compared 
to the lumped parameter analysis for free-free 
boundaries in Figures 9,10 and 11 for various 
models. 


AMPLITUDE AMPLITUDE 


propagation nathod 

ltjupad-paraaiatar nathod 



0 100 200 


FRBQUSKOY (kH*) 

Fig. 9 - Composite cone, natural frequencies, 
polyethylene base, I uc i te tip 


war* propagation method 

lunpod-paramater nathod 



FREQUENCY (kE») 

Fig. II - Composite cone, natural frequen- 
cies, polyethylene shell base, 
polyethylene tip 

DISCUSSION OF RESULTS 

The natural frequencies for various com- 
posite rods, cones and conical shells obtained 
experimentally and by the lumped parameter 
approach compare favorably. The maximum de- 
viation is less than 3% for the composite cones 
and conical shells. The second natural fre- 
quency for the pol ye thy lene-1 uc i te cone was 
not obtained experimentally because of 
difficulties in the experimental set-up(Table 
3). Subsequent modifications eliminated the 
difficulties, however the dzta was not retaken. 


prop*s«tion natXod 

lunpod-porasater uthod 



FREQUENCY (kH*) 

Fig. 10 - Composite cone, natural frequen- 
-cies, lucite shell base, 
polyethylene tip 


Unfortunately the material properties of 
the different plastic materials do not vary 
greatly. The use of materials having a wider 
range of properties should give more sig- 
nificant results for composite bodies. Metals 
offer certain advantages over plastics such as 
a wider range of properties, more stable proper- 
ties, and easier machi ni ng-part icul ar 1 y for 
conical shells. The use of metals is recom- 
mended if a higher frequency input is available. 

The mode shapes and relative strains are 
presented from the lumped parameter approach 
only. Since strain gages were employed at only 
two locations on each model, comparison was 
possible only at those points. Therefore these 
results were not considered significant and are 
not presented. 

Comparison of the two theoretical pre- 
dictions, lumped parameter and wave propagation 
for the free-free boundary was very erratic and 
poor. The fundamental natural frequencies for 
the various composite cones and conical shells 


233 



compared reasonably well but not the higher 
modes. This confirmed the fact that while 
elastic wave propagation analysis might give 
good results for problems with a short duration 
input loading, it does not give very good 
results for long time input loading. In 
addition, the wave propagation analysis is a 
one-dimensional analysis. For conical modes 
where the area varies, it is dependent on the 
cone angle and the number of grid points 
used in representing the model. For one 
dimensional problems such as longitudinal 
vibration of composite rods, axially 
symmetric longitudinal vibration of circular 
cylinders and spherical shells, good agree- 
.ments^can. be expected. (Natural frequencies 

were obtained for two-layer composite rods, 
cylindrical and spherical shells using the wave 
propagation approach. Results compared well 
with the theoretically computed results. The 
results for these calculations and comparison 
are not presented in this paper.) 

CONCLUSIONS 

< 

The lumped parameter method using a series 
of mass elements and linear springs can be used 
to analyze the longitudinal vibration problem 
of composite ax i symmetric bodies of varying 
area as verified by experimental vibration 
tests on simple composite rods, cones and con- 
ical shells. 

The one dimensional elastic wave prop- 
agation approach utilized in computing the 
natural frequencies of the composite models 
did not give good results. In order to obtain 
better results using the one dimensional wave 
propagation analysis for the changing area 
problem with a long time input loading, much 
finer gird points must be utilized which re- 
quire increased running time and storage space 
on the computer. 

More research is needed to determine the 
effect of cone angle, impedance change, mater- 
ial properties , etc. on the vibration charac- 
teristic. 
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